# Replication package for 
# "International Agreement Design and the Moderating Role of Domestic Bureaucratic Quality: The Case of Freshwater Cooperation"
# Johannes Karreth and Jaroslav Tir
# Forthcoming in the Journal of Peace Research
# November 2017
# jkarreth@ursinus.edu

# Please read "0_readme.html" before using this code.

######################################
### Load packages
######################################

rm(list = ls())

library("foreign")
library("texreg")
library("lme4")
library("xtable")
library("dplyr")
library("gridExtra")
library("effects")
library("ggplot2")
library("rio")

source("ggintfun.R")

# Set working directory
# setwd("...")
setwd("/Users/johanneskarreth/Documents/Dropbox/Uni/1 - Papers/34 - River treaties/Paper/Project/JPR submission/Final/Replication")

######################################
### Load data
######################################

dat <- import("basins_small.csv")

dat <- dat[, c("year", "dyadid", "ccode1", "ccode2", "baravg2", "numbasins", "instcoop_max", "instcoop_min", "instcoop_median", "instcoop_mean", "ntreaties", "ntreaties0", "ntreaties2", "instcoop_mean2", "instcoop_max2", "contiguity", "lnpowerratio", "interdepGl", "interdepBa", "dyaddem6", "dyaddem7", "icrg_qog_min", "icrg_qog_max", "icrg_burqual_min", "icrg_demacc_min", "alliance", "depratio_max", "depratio_min", "wres_min_log", "totwithdraw_min_log", "gdpmax_log", "gdpmaxPWT_log", "monitoring_max", "enforcement_max", "conflictres_max", "institution_max", "rivalry_th", "cie_min", "cie_bin_min", "instcoop_max2_lag1", "icrg_burqual_min_lag1", "icrg_qog_min_lag1", "numstate")]

######################################
### Add BAR data
######################################
  
bar <- import("EventMaster111710.csv")

bar <- bar[, c("DATE", "DYAD_CODE", "BCode", "NUMBER_OF_BASINS", "LOCATION", "CCODE1", "CCODE2", "BCCODE1", "BCCODE2", "NUMBER_OF_Countries", "COPDAB_SCALE", "BAR_Scale", "Issue_Type1", "Issue_Type2", "EVENT_ISSUE")]

library("countrycode")
bar$ccode1 <- countrycode(bar$CCODE1, origin = "iso3c", destination = "cown")
bar$ccode2 <- countrycode(bar$CCODE2, origin = "iso3c", destination = "cown")

# Add country codes for countries not properly matched

bar[bar$CCODE1 == "CZS", ]$ccode1 <- 315
bar[bar$CCODE1 == "DRV", ]$ccode1 <- 816
bar[bar$CCODE1 == "GDR", ]$ccode1 <- 265  
bar[bar$CCODE1 == "GFR", ]$ccode1 <- 260  
bar[bar$CCODE1 == "PRK", ]$ccode1 <- 731
bar[bar$CCODE1 == "ROM", ]$ccode1 <- 360
bar[bar$CCODE1 == "SRB", ]$ccode1 <- 345
bar[bar$CCODE1 == "USR", ]$ccode1 <- 365
bar[bar$CCODE1 == "YGF", ]$ccode1 <- 345  
bar[bar$CCODE1 == "ZAR", ]$ccode1 <- 490  

bar[bar$CCODE2 == "CZS", ]$ccode2 <- 315
bar[bar$CCODE2 == "DRV", ]$ccode2 <- 816
bar[bar$CCODE2 == "GDR", ]$ccode2 <- 265  
bar[bar$CCODE2 == "GFR", ]$ccode2 <- 260
bar[bar$CCODE2 == "PRK", ]$ccode2 <- 731
bar[bar$CCODE2 == "ROM", ]$ccode2 <- 360
bar[bar$CCODE2 == "SRB", ]$ccode2 <- 345
bar[bar$CCODE2 == "USR", ]$ccode2 <- 365
bar[bar$CCODE2 == "YGF", ]$ccode2 <- 345  
bar[bar$CCODE2 == "ZAR", ]$ccode2 <- 490  
  
bar$CCODE1 <- NULL
bar$CCODE2 <- NULL

bar <- bar[!is.na(bar$ccode1) == TRUE, ]
bar <- bar[!is.na(bar$ccode2) == TRUE, ]

bar$dyadid <- with(bar, ifelse(ccode1 < ccode2, ccode1 * 1000 + ccode2, ccode2 * 1000 + ccode1))

bar$date <- as.Date(bar$DATE, format = "%m/%d/%y", origin = "1900-01-01")
bar$year <- as.numeric(format(bar$date, "%Y"))
bar$year <- ifelse(bar$year > 2008, (bar$year - 100), bar$year)

bar <- bar[!is.na(bar$year) == TRUE, ]

bar$DATE <- NULL
bar$date <- NULL

bar$bar2 <- with(bar, ifelse(is.na(BAR_Scale) == TRUE, 0, BAR_Scale))
bar$bar2_pos <- with(bar, ifelse(bar2 > 0, 1, 0))
bar$bar2_neg <- with(bar, ifelse(bar2 < 0, 1, 0))
bar.dy <- summarise(group_by(bar, dyadid, year, add = FALSE), 
                    bar_pos = sum(bar2_pos), 
                    bar_neg = sum(bar2_neg), 
                    bar_avg = mean(bar2), 
                    bar_median = quantile(bar2, probs = c(0.5)),
                    bar_pos_avg = mean(bar2_pos),
                    bar_neg_avg = mean(bar2_neg))
bar.dy <- arrange(bar.dy, dyadid, year)

dat <- merge(x = dat, y = bar.dy, by = c("dyadid", "year"), all.x = TRUE, all.y = FALSE)
  
######################################
### Cleaning
######################################

dat <- dat[dat$year > 1983, ]
summary(dat)

### _2 variables fill in 0s for all years with no events, like so:
### bar_pos: only for years with events, non-event years missing
### bar_pos2: for all years, non-event years coded as 0

dat$bar_pos2 <- ifelse(is.na(dat$bar_pos) == TRUE, 0, dat$bar_pos)
dat$bar_neg2 <- ifelse(is.na(dat$bar_neg) == TRUE, 0, dat$bar_neg)
dat$bar_pos_avg2 <- ifelse(is.na(dat$bar_pos_avg) == TRUE, 0, dat$bar_pos_avg)
dat$bar_neg_avg2 <- ifelse(is.na(dat$bar_neg_avg) == TRUE, 0, dat$bar_neg_avg)
dat$bar_pos2only <- ifelse(dat$bar_pos2 > 0 & dat$bar_neg2 == 0, 1, 0) # Also need version for years w/ events only
dat$bar_pos2only2 <- ifelse(is.na(dat$bar_avg) == FALSE, ifelse(dat$bar_pos2 > 0 & dat$bar_neg2 == 0, 1, 0), NA)
dat$bar_neg2only <- ifelse(dat$bar_neg2 > 0 & dat$bar_pos2 == 0, 1, 0) # Also need version for years w/ events only
dat$bar_neg2only2 <- ifelse(is.na(dat$bar_avg) == FALSE, ifelse(dat$bar_neg2 > 0 & dat$bar_pos2 == 0, 1, 0), NA)
dat$bar_neg2any <- ifelse(dat$bar_neg2 > 0 & is.na(dat$bar_neg2) == FALSE, 1, dat$bar_neg2)
dat$bar_avg2 <- ifelse(is.na(dat$bar_avg) == TRUE, 0, dat$bar_avg)
dat <- data.frame(mutate(group_by(dat, dyadid), bar_avg2.l1 = lag(bar_avg2, n = 1)))
dat$bar_median2 <- ifelse(is.na(dat$bar_median) == TRUE, 0, dat$bar_median)

# Rescale the bureaucratic quality measure to 0-1
range01 <- function(x, na.rm = TRUE){(x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))}
dat$icrg_burqual_min_std <- range01(dat$icrg_burqual_min, na.rm = TRUE)

### Alternative Institutionalization coding: instcoop_max_05
## 0=no treaty, 1=treaty w/o institutional features, 2=1 institutional feature, etc. Result: 0-5 scale

dat$instcoop_max_05 <- dat$instcoop_max + 1
dat[dat$ntreaties0 == 1, ]$instcoop_max_05 <- 0
table(dat[dat$ntreaties0 == 0,]$instcoop_max2)
table(dat$instcoop_max_05, dat$ntreaties0)

### 4 treaty components: create alternative versions where no info = 0
## Similar to instcoop_max2

dat$institution_max2 <- ifelse(is.na(dat$institution_max), 0, dat$institution_max)
dat$monitoring_max2 <- ifelse(is.na(dat$monitoring_max), 0, dat$monitoring_max)
dat$enforcement_max2 <- ifelse(is.na(dat$enforcement_max), 0, dat$enforcement_max)
dat$conflictres_max2 <- ifelse(is.na(dat$conflictres_max), 0, dat$conflictres_max)

######################################
### Describe data
######################################

### Spread of treaty institutionalization into LDCs

rt.dat <- import("rtd_treaty-country-year.csv")
rt.dat <- rt.dat[, c("cowcode", "year", "treaty_metadata_docid", "treatyinst.ind", "treatyinst.ind2")]
rt.dat$treaty <- c(rep(1, length(rt.dat$cowcode)))
rt.dat.country <- summarize(group_by(rt.dat, cowcode, year), treatycount = sum(treaty), treatyinst.ind.max = max(treatyinst.ind), treatyinst.ind2.max = max(treatyinst.ind2), treatyinst.ind.min = min(treatyinst.ind), treatyinst.ind2.min = min(treatyinst.ind2))

oecd <- import("oecd_members.csv")
rt.dat.country <- merge(x = rt.dat.country, y = oecd, by.x = c("cowcode"), by.y = c("cowcode"), all.x = TRUE, all.y = FALSE)
rt.dat.country$oecd <- ifelse(rt.dat.country$year >= rt.dat.country$startyear, 1, 0)
rt.dat.country$oecd <- ifelse(is.na(rt.dat.country$oecd) == TRUE, 0, rt.dat.country$oecd)

data.frame(summarize(group_by(rt.dat.country[rt.dat.country$year >= 1961, ], year), 
          treatyyes_nonOECD = 1 - mean(oecd)))

rt.dat.country$region <- countrycode::countrycode(rt.dat.country$cowcode, origin = "cown", destination = "region")
rt.dat.country.region <- summarize(group_by(rt.dat.country, region, year), 
                     avg_treatycount = mean(treatycount),
                     avg_treatyinst = mean(treatyinst.ind.max))
rt.dat.country.region[rt.dat.country.region$year == 1960, ]
rt.dat.country.region[rt.dat.country.region$year == 2000, ]

library("directlabels")
p1 <- ggplot(data = rt.dat.country.region[rt.dat.country.region$region != "Melanesia" & rt.dat.country.region$year < 2007, ], aes(x = year, y = avg_treatycount, group = region, color = region)) + geom_line() + theme_jk() + xlab("") + ylab("Average treaties per dyad") + xlim(1945,2010) + geom_vline(xintercept = 1984) + geom_vline(xintercept = 2006) 
p2 <- direct.label(p1)
ggsave(p2, file = "Output/rt_count_spread.pdf", width = 7, height = 5)

rt.dat.country$continent <- countrycode::countrycode(rt.dat.country$cowcode, origin = "cown", destination = "continent")
rt.dat.country.continent <- summarize(group_by(rt.dat.country, continent, year), 
                                      avg_treatycount = mean(treatycount),
                                      avg_treatyinst = mean(treatyinst.ind.max))
rt.dat.country.continent[rt.dat.country.continent$year == 1960, ]
rt.dat.country.continent[rt.dat.country.continent$year == 2000, ]

p3 <- ggplot(data = rt.dat.country.continent[rt.dat.country.continent$continent != "Oceania" & is.na(rt.dat.country.continent$continent) == FALSE & rt.dat.country.continent$year > 1960 & rt.dat.country.continent$year < 2007, ], aes(x = year, y = avg_treatyinst, group = continent, color = continent)) + geom_line() + theme_jk() + xlab("") + ylab("Average treaty institutionalization per dyad") + xlim(1960, 2010) + geom_vline(xintercept = 1984) + geom_vline(xintercept = 2006)
p4 <- direct.label(p3, method = "last.points")

ggsave(p4, file = "Output/rt_inst_spread.pdf", width = 7, height = 5)

### Bar event codes

codes <- data.frame(table(as.integer(dat$bar_avg)))
codes$desc <- c("Small-scale military acts", "Political/military hostile acts", "Diplomatic/economic hostile acts", "Strong/official verbal hostility", "Mild/unofficial verbal hostility", "Neutral, non-significant acts", "Mild verbal support", "Official verbal support", "Cultural, scientific agreement/support", "Non-military economic, technological, industrial agreement", "Military, economic, strategic support", "International water treaty")
colnames(codes) <- c("Code", "Frequency", "Description")
codes$Code <- round(as.numeric(levels(codes$Code)), 0)

codes.tab <- xtable(codes[, c(1, 3)], caption = "", label = "tab:BAR_codes", align = "lll", digits = 0)
# print(codes.tab, file = "Output/barcodes.tex", table.placement = "hbtp", caption.placement = "top", size = "footnotesize", include.rownames = FALSE, include.colnames = TRUE, math.style.negative = TRUE, booktabs = TRUE)

### List of key variables for discussion of case evidence

dat.cases <- dat[, c("dyadid", "ccode1", "ccode2", "year", "bar_avg2", "instcoop_max2", "icrg_burqual_min_std", "wres_min_log", "ntreaties2", "dyaddem7", "gdpmaxPWT_log", "lnpowerratio", "alliance")]
dat.cases <- summarize(group_by(dat.cases, dyadid),
                       ccode1 = ccode1[1],
                       ccode2 = ccode2[2],
                       Avg_Cooperation = mean(bar_avg2, na.rm = TRUE), 
                       Avg_Treaty_Institutionalization = mean(instcoop_max2, na.rm = TRUE), 
                       Avg_Lower_BQ = mean(icrg_burqual_min_std, na.rm = TRUE), 
                       Avg_Water_Availability = mean(wres_min_log, na.rm = TRUE), 
                       Avg_Number_Treaties = mean(ntreaties2, na.rm = TRUE), 
                       Dem_Dyad = mean(dyaddem7, na.rm = TRUE), 
                       Avg_GDPpc_max = mean(gdpmaxPWT_log, na.rm = TRUE), 
                       Avg_Power_Ratio = mean(lnpowerratio, na.rm = TRUE), 
                       Alliance = mean(alliance, na.rm = TRUE))

dat.cases$country1 <- countrycode::countrycode(dat.cases$ccode1, origin = "cown", destination = "country.name")
dat.cases$country2 <- countrycode::countrycode(dat.cases$ccode2, origin = "cown", destination = "country.name")

dat.cases <- dat.cases[, c(13, 14, 4:12, 1, 2, 3)]
dat.cases <- dat.cases[is.na(dat.cases$Avg_Lower_BQ) == FALSE & is.na(dat.cases$Avg_Treaty_Institutionalization) == FALSE, ]

export(dat.cases, file = "Output/RT_Dyads_Averages.xlsx")

dat.cases <- dat[, c("dyadid", "ccode1", "ccode2", "year", "bar_avg2", "instcoop_max2", "icrg_burqual_min_std", "wres_min_log", "ntreaties2", "dyaddem7", "gdpmaxPWT_log", "lnpowerratio", "alliance")]
dat.cases <- summarize(group_by(dat.cases[dat.cases$year > 2000, ], dyadid),
                       ccode1 = ccode1[1],
                       ccode2 = ccode2[2],
                       Avg_Cooperation = mean(bar_avg2, na.rm = TRUE), 
                       Avg_Treaty_Institutionalization = mean(instcoop_max2, na.rm = TRUE), 
                       Avg_Lower_BQ = mean(icrg_burqual_min_std, na.rm = TRUE), 
                       Avg_Water_Availability = mean(wres_min_log, na.rm = TRUE), 
                       Avg_Number_Treaties = mean(ntreaties2, na.rm = TRUE), 
                       Dem_Dyad = mean(dyaddem7, na.rm = TRUE), 
                       Avg_GDPpc_max = mean(gdpmaxPWT_log, na.rm = TRUE), 
                       Avg_Power_Ratio = mean(lnpowerratio, na.rm = TRUE), 
                       Alliance = mean(alliance, na.rm = TRUE))

dat.cases$country1 <- countrycode::countrycode(dat.cases$ccode1, origin = "cown", destination = "country.name")
dat.cases$country2 <- countrycode::countrycode(dat.cases$ccode2, origin = "cown", destination = "country.name")

dat.cases <- dat.cases[, c(13, 14, 4:12, 1, 2, 3)]
dat.cases <- dat.cases[is.na(dat.cases$Avg_Lower_BQ) == FALSE & is.na(dat.cases$Avg_Treaty_Institutionalization) == FALSE, ]

export(dat.cases, file = "Output/RT_Dyads_post2000.xlsx")

### Aggregation: positive counts (versus years with no events)

desc.dat <- na.omit(dat[, c("bar_pos2only", "bar_avg2", "bar_pos2", "instcoop_max2", "instcoop_max_05", "rivalry_th", "icrg_qog_min", "icrg_burqual_min", "icrg_burqual_min_std", "wres_min_log", "ntreaties2", "dyaddem7", "gdpmaxPWT_log", "lnpowerratio", "alliance", "dyadid", "year")])

# Counts of positive eventes
table(desc.dat[, "bar_pos2"])

# Positive events only
table(desc.dat[, "bar_pos2only"])

### Aggregation: yearly averages

table(as.integer(desc.dat[, "bar_avg2"]))
table(cut(x = desc.dat[, "bar_avg2"], breaks = c(min(desc.dat[, "bar_avg2"]), -0.0001, 0.0001, max(desc.dat[, "bar_avg2"])), include.lowest = TRUE))

### Variance across years and dyads

# By year
barmeans.year <- summarise(group_by(desc.dat, year), bar_avg = mean(bar_avg2), bar_sd = sd(bar_avg2) / sqrt(length(bar_avg2)))

p.barmeans.year <- ggplot(data = barmeans.year, aes(x = year, y = bar_avg)) + geom_line(color = "gray")  + geom_point() + geom_segment(aes(x = year, xend = year, y = bar_avg - bar_sd, yend = bar_avg + bar_sd))+ theme_jk() + xlab("Year") + ylab("Average conflict/cooperation score across all dyads") + scale_x_continuous(limits = c(1984, 2006), breaks = c(1984, 1990, 1995, 2000, 2006), labels = c(1984, 1990, 1995, 2000, 2006))

 # By dyad - all years
barmeans.dyadid <- summarise(group_by(desc.dat, dyadid), bar_avg = mean(bar_avg2), bar_sd = sd(bar_avg2) / sqrt(length(bar_avg2)), obs = length(bar_avg2))
barmeans.dyadid$dyadid2 <- as.numeric(as.ordered(barmeans.dyadid$dyadid))

p.barmeans.dyadid <- ggplot(data = barmeans.dyadid, aes(x = reorder(dyadid2, bar_avg), y = bar_avg)) + geom_segment(aes(x = reorder(dyadid2, bar_avg), xend = reorder(dyadid2, bar_avg), y = bar_avg - bar_sd, yend = bar_avg + bar_sd, color = obs), alpha = 0.5) + geom_point(aes(color = obs)) + ylim(c(-1, 3)) + guides(alpha = FALSE, color = FALSE) + coord_flip() + scale_x_discrete(breaks = NULL, expand = c(0.01, 0.01)) + theme_jk() + xlab("Dyads") + ylab("Average conflict/cooperation score across all dyad-years")
p.barmeans.dyadid <- ggplot(data = barmeans.dyadid, aes(x = dyadid2, y = bar_avg)) + geom_segment(aes(x = dyadid2, xend = dyadid2, y = bar_avg - bar_sd, yend = bar_avg + bar_sd, color = obs), alpha = 0) + geom_point(aes(color = obs)) + ylim(c(-0.25, 2.25)) + guides(alpha = FALSE, color = FALSE) + coord_flip() + scale_x_discrete(breaks = NULL, expand = c(0.01, 0.01)) + theme_jk() + xlab("Dyads") + ylab("Average conflict/cooperation score across all dyad-years")

p.barmeans.dyadid <- ggplot(data = barmeans.dyadid, aes(x = dyadid2, y = bar_avg)) + geom_segment(aes(x = dyadid2, xend = dyadid2, y = bar_avg - bar_sd, yend = bar_avg + bar_sd), alpha = 0.5) + geom_point(size = 1) + ylim(c(-1, 3)) + guides(alpha = FALSE, color = FALSE) + scale_x_discrete(breaks = NULL, expand = c(0.01, 0.01)) + theme_jk() + xlab("Dyads") + ylab("Average conflict/cooperation score across all dyad-years")

pdf("Output/barmeans.pdf", width = 12, height = 6)
grid.arrange(p.barmeans.year + ggtitle("Variance across time"), p.barmeans.dyadid + ggtitle("Variance across dyads"), ncol = 2)
dev.off()

### Histograms of treaty institutionalization

desc.dat$icrg_burqual_min_q <- as.factor(cut(x = desc.dat$icrg_burqual_min, breaks = quantile(desc.dat$icrg_burqual_min, probs = c(0.1, 0.5, 0.9, 1), na.rm = TRUE), include.lowest = TRUE))
desc.dat$icrg_burqual_min_bin10 <- as.numeric(cut(desc.dat$icrg_burqual_min_std, breaks = 10, include.lowest = TRUE, right = TRUE))

levels(desc.dat$icrg_burqual_min_q) <- c("BQ: 10th percentile", "BQ: Median", "BQ: 90th percentile")

inst.hist <- ggplot(data = desc.dat[desc.dat$instcoop_max2 > 0, ], aes(x = as.factor(instcoop_max2), group = icrg_burqual_min_q, fill = icrg_burqual_min_q)) + geom_bar(stat = "count", position = "dodge", alpha = 0.5) + theme_jk() + xlab("Treaty institutionalization") + ylab("Dyad-years") + facet_wrap(~ icrg_burqual_min_q, scales = "free_y") + guides(fill = FALSE)

pdf("Output/insthist.pdf", width = 6, height = 3)
inst.hist
dev.off()

### Histograms of Bureaucratic quality

bq.hist <- ggplot(data = data.frame(table(cut(desc.dat$icrg_burqual_min_std, breaks = 5))), aes(x = as.numeric(Var1), y = Freq)) + geom_bar(stat = "identity", width = 1, fill = "darkgray", alpha = 0.75) + scale_x_continuous(breaks = seq(0.5, 5.5, by = 1), labels = seq(0, 1, by = 0.2)) + theme_jk() + xlab("Bureaucratic quality") + ylab("Dyad-years")

pdf("Output/bqhist.pdf", width = 4, height = 4)
bq.hist
dev.off()


######################################
### Diagnostics for common support (not in paper)
######################################

# All cases

dat$icrg_burqual_min_tri <- cut(dat$icrg_burqual_min_std, breaks = c(0, 0.25, 0.5, 1), include.lowest = TRUE)
dat$icrg_burqual_min_tri <- ifelse(dat$icrg_burqual_min_tri == "(0.25,0.5]", "BQ = [0.25, 0.5]",
                                   ifelse(dat$icrg_burqual_min_tri == "(0.5,1]", "BQ = [0.5, 1]",
                                          ifelse(dat$icrg_burqual_min_tri == "[0,0.25]", "BQ = [0, 0.25]", NA)))
diag_p_avg2 <- ggplot(data = dat[is.na(dat$icrg_burqual_min_tri) == FALSE, ], aes(x = jitter(instcoop_max2), y = jitter(bar_avg2))) + geom_point(shape = 20, size = 0.5) + geom_smooth(method = "loess", color = "red", size = 0.5, se = FALSE) + geom_smooth(method = "lm", color = "blue", size = 0.5, se = FALSE) + facet_wrap(~ icrg_burqual_min_tri, ncol = 3) + xlab("Treaty institutionalization") + ylab("Observed cooperation") + theme_jk()

# Years with events only

dat$icrg_burqual_min_tri2 <- ifelse(is.na(dat$bar_avg) == FALSE, cut(dat$icrg_burqual_min_std, breaks = c(0, 0.25, 0.5, 1), include.lowest = TRUE), NA)
dat$icrg_burqual_min_tri2 <- ifelse(dat$icrg_burqual_min_tri2 == 2, "BQ = [0.25, 0.5]",
                                   ifelse(dat$icrg_burqual_min_tri2 == 3, "BQ = [0.5, 1]",
                                          ifelse(dat$icrg_burqual_min_tri2 == 1, "BQ = [0, 0.25]", NA)))
diag_p_avg <- ggplot(data = dat[is.na(dat$icrg_burqual_min_tri2) == FALSE, ], aes(x = jitter(instcoop_max2), y = jitter(bar_avg))) + geom_point(shape = 20, size = 0.5) + geom_smooth(method = "loess", color = "red", size = 0.5, se = FALSE) + geom_smooth(method = "lm", color = "blue", size = 0.5, se = FALSE) + facet_wrap(~ icrg_burqual_min_tri2, ncol = 3) + xlab("Treaty institutionalization") + ylab("Observed cooperation") + theme_jk()

pdf("Output/ia_diag_both.pdf", width = 9, height = 8)
grid.arrange(diag_p_avg2 + ggtitle("All years"), diag_p_avg + ggtitle("Years with events only"), ncol = 1)
dev.off()


######################################
### Analyses
######################################

### Main analyses

# 1.1 OLS
ols.bq <- lm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

# 1.2 Probit
probit.bq <- glm(bar_pos2only ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance , data = dat, x = TRUE, family = binomial(link = "probit"))


### Robustness tests

# 2. Years with events only

ols.bq.eo <- lm(bar_avg ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

# 3. Median cooperation instead of mean/average

ols.bq.median <- lm(bar_median2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

# 4. Indicator for no treaty

# Add-on: 3-way interaction with treaty indicator
dat$notreaties <- ifelse(dat$ntreaties0 == 1, 0, 1)

ols.bq.3w <- lm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std * notreaties + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

# 5. Dyad FE

ols.bq.dfe <- lm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance + factor(dyadid), data = dat, x = TRUE)

# 6a. Dyad FELDV
# Note: I'm aware that the estimates below are likely biased (because FEs are correlated with bar_avg2.l1),
# but this robustness test follows a reviewer's request.
# As an alternative, I also show results from using the Arellano-Bond estimator below.

ols.bq.dfeldv <- lm(bar_avg2 ~ bar_avg2.l1 + instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance + factor(dyadid), data = dat, x = TRUE)

# 6b. Arellano-Bond GMM

library("plm")
pdat <- pdata.frame(na.omit(dat[, c("bar_avg2", "instcoop_max2", "icrg_burqual_min_std", "wres_min_log", "ntreaties2", "dyaddem7", "gdpmaxPWT_log", "lnpowerratio", "alliance", "dyadid", "year")]), index = c("dyadid", "year"), drop.index = FALSE, row.names = TRUE)
ols.abgmm <- pgmm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance | lag(bar_avg2, 2:99), data = pdat, effect = "individual")
detach("package:plm")

# 7. Year FE

ols.bq.yfe <- lm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance + factor(year), data = dat, x = TRUE)

# 8. Dyad & year RE

ols.bq.mlm <- lmer(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance + (1 | dyadid) + (1 | year), data = dat)

# 9. Rivalry

ols.bq.riv <- lm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance + rivalry_th, data = dat)

# 10. IV

export(dat[, c("bar_avg2", "instcoop_max2", "icrg_burqual_min_std", "wres_min_log", "ntreaties2", "dyaddem7", "gdpmaxPWT_log", "lnpowerratio", "alliance", "numstate", "dyadid", "year")], file = "Output/rt_iv_data.dta", version = 13)

# Use the RStata package to access Stata directly (or execute the code in Stata)
library("RStata")
options("RStata.StataPath" = "/Applications/Stata/StataSE.app/Contents/MacOS/stata-se")
options("RStata.StataVersion" = 13)

stata_cmd <- "
use rt_iv_data.dta
quietly reg bar_avg2  c.instcoop_max2##c.icrg_burqual_min_std wres_min_log ntreaties2 dyaddem7 gdpmaxPWT_log lnpowerratio alliance
gen burqual_sample = 1 if e(sample)
set matsize 5000
quietly reg instcoop_max2 numstate i.dyadid i.year if burqual_sample == 1
predict instcoop_max2_IV, xb
gen instcoopXburqual = instcoop_max2 * icrg_burqual_min
gen instcoopXburqual_IV = instcoop_max2_IV * icrg_burqual_min
ivregress 2sls bar_avg2 icrg_burqual_min_std wres_min_log ntreaties2 dyaddem7 gdpmaxPWT_log lnpowerratio alliance (instcoop_max2 instcoopXburqual = instcoop_max2_IV instcoopXburqual_IV)
matrix list e(V)
"
stata(stata_cmd)   

# 11. Breaking up institutions

ols.bq.institution <- lm(bar_avg2 ~ institution_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

ols.bq.monitoring <- lm(bar_avg2 ~ monitoring_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

ols.bq.enforcement <- lm(bar_avg2 ~ enforcement_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

ols.bq.conflictres <- lm(bar_avg2 ~ conflictres_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

# 12. Seemingly unrelated regression

library("systemfit")

eq.baravg <- bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance
eq.instcoop <- instcoop_max2 ~ icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance


ols.baravg <- lm(eq.baravg, data = dat)
ols.instcoop <- lm(eq.instcoop, data = dat)

sur.bq <- systemfit(list(lm.instcoop = eq.instcoop, lm.baravg = eq.baravg), data = dat)

# Also show plot:

inst_bq_mod <- lm(instcoop_max2 ~ icrg_burqual_min_std, data = dat)

inst_bq_p <- ggplot(data = dat, aes(x = icrg_burqual_min_std, y = instcoop_max2)) + geom_jitter(width = 0.15, height = 0.5, shape = 3, size = 0.25, color = "grey") + geom_smooth(method = "loess", fullrange = TRUE, color = "black", fill = "gray", linetype = "dashed", size = 0.75, se = FALSE) + geom_smooth(method = "lm", fullrange = TRUE, color = "black", fill = "gray", size = 0.75, se = FALSE) + theme_jk() + scale_x_continuous(breaks = c(seq(0, 1, by = 0.2)), labels = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + xlab("Bureaucratic quality") + ylab("Treaty institutionalization") + ggtitle(paste("R-squared = ", signif(summary(inst_bq_mod)$r.squared, 1), ", Slope = ",signif(inst_bq_mod$coef[[2]], 1), ", P < 0.001", sep = ""))

ggsave(inst_bq_p, filename = "Output/inst_bq.pdf", width = 5, height = 4)

# 13. Control for shared IGO memberships

igo.dat <- import("IGO_dyadunit_v2.3.RData")
igo.dat$sharedIGOcount <- rowSums(igo.dat[, c(6:533)] == 1, na.rm = TRUE)
igo.dat <- igo.dat[igo.dat$year > 1983, c("ccode1", "ccode2", "year", "sharedIGOcount")]
igo.dat$dyadid <- ifelse(igo.dat$ccode1 < igo.dat$ccode2, igo.dat$ccode1 * 1000 + igo.dat$ccode2, igo.dat$ccode2 * 1000 + igo.dat$ccode1)
dat <- merge(x = dat, y = igo.dat, by.x = c("dyadid", "year"), by.y = c("dyadid", "year"), all.x = TRUE, all.y = FALSE)

ols.bq.igos <- lm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance + sharedIGOcount, data = dat, x = TRUE)

# 14. Conflict events

ols.bq.neg2 <- lm(bar_neg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

ols.bq.negavg2 <- lm(bar_neg_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE)

probit.bq.neg2any <- glm(bar_neg2any ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance, data = dat, x = TRUE, family = binomial(link = "probit"))

mid.dat <- import("dyadicmid2.0.dta")
mid.dat$mid <- 1
mid.dat$dyadid <- ifelse(mid.dat$statea < mid.dat$stateb, mid.dat$statea * 1000 + mid.dat$stateb, mid.dat$stateb * 1000 + mid.dat$statea)
dat.mid <- merge(x = dat, y = mid.dat[, c("dyadid", "year", "mid", "fatlev")], by.x = c("dyadid", "year"), by.y = c("dyadid", "year"), all.x = TRUE, all.y = FALSE)
dat.mid$mid <- ifelse(is.na(dat.mid$mid) == TRUE & dat.mid$year < 2002, 0, dat.mid$mid)

probit.bq.mid <- glm(mid ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance , data = dat.mid, x = TRUE, family = binomial(link = "probit"))

# 15. States in treaty

rtd <- import("rtd_dyad-year.csv")
dat <- merge(x = dat, y = rtd, by = c("dyadid", "year"), all.x = TRUE, all.y = FALSE)
dat$treatymembers.max2 <- ifelse(is.na(dat$treatymembers.max) == TRUE, 0, dat$treatymembers.max)
ols.bq.tm <- lm(bar_avg2 ~ instcoop_max2 * icrg_burqual_min_std + wres_min_log + ntreaties2 + dyaddem7 + gdpmaxPWT_log + lnpowerratio + alliance + treatymembers.max2, data = dat, x = TRUE)

######################################
### Tables
######################################

### Table 1 (main results)

tab1 <- texreg(list(ols.bq, probit.bq), 
               file = "Output/tab_main1.tex", 
               stars = 0.05, 
               custom.model.names = c("OLS", "Probit"), 
               custom.coef.names = c("Intercept", "Treaty institutionalization", "Bureaucratic quality (lower)", "Water availability (lower)", "Treaty count", "Democratic dyad", "GDP p.c. (higher, logged)", "Power ratio", "Alliance", "Treaty institutionalization $times$ Bureaucratic quality (lower)"), 
               custom.note = "%stars. Standard errors in parentheses.", 
               digits = 3, 
               leading.zero = TRUE, 
               # omit.coef = "dyadid", 
               reorder.coef = c(2, 3, 10, 4, 5, 6, 7, 8, 9, 1),
               center = TRUE, 
               caption = "Determinants of water-related cooperation.", 
               caption.above = TRUE, 
               label = "tab:main", 
               dcolumn = TRUE, 
               booktabs = TRUE, 
               use.packages = FALSE, 
               table = TRUE)
               

### Table 2 (first set of robustness tests)

# Use ols.bq as placeholder for IV results
# then replace later by hand
# (don't forget to update SEs & p-value markers, too!!!)

tab2 <- texreg(list(ols.bq.median, ols.bq.eo, ols.bq.3w, ols.bq.dfe, ols.bq.dfeldv, ols.abgmm, ols.bq.yfe, ols.bq.mlm, ols.bq.riv, ols.bq, ols.bq.igos, ols.bq.tm),
		  	   file = "Output/tab_robust1.tex", 
               stars = 0.05, 
               # custom.model.names = c("1", "2", "1-EO"), 
               # custom.coef.names = c("Intercept", "Treaty institutionalization", "Bureaucratic quality (lower)", "Water availability (lower)", "Treaty count", "Democratic dyad", "GDP p.c. (higher, logged)", "Power ratio", "Alliance", "Treaty institutionalization $times$ Bureaucratic quality (lower)", "No treaties", "No treaties $times$ Bureaucratic quality (lower)", "Cooperation (previous year)", "Rivalry", rep("foo", 761)), 
               custom.note = "%stars. Standard errors in parentheses.", 
               digits = 3, 
               leading.zero = TRUE, 
               omit.coef = "factor", 
               reorder.coef = c(2, 3, 10, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 15, 16, 1),
               center = TRUE, 
               caption = "Determinants of water-related cooperation: Robustness tests. UPDATE IV ESTIMATES BY HAND!!!", 
               caption.above = TRUE, 
               label = "tab:robust1", 
               dcolumn = TRUE, 
               booktabs = TRUE, 
               use.packages = FALSE, 
               table = TRUE)

### Table 3 (breaking up the institutions index)
# "Treaty created IGO", "Treaty provides Monitoring", "Treaty provides Enforcement", "Treaty provides conflict resolution"

tab3 <- texreg(list(ols.bq.institution, ols.bq.monitoring, ols.bq.enforcement, ols.bq.conflictres),
		  	   file = "Output/tab_robust2.tex", 
               stars = 0.05, 
               # custom.model.names = c("1", "2", "1-EO"), 
               custom.coef.names = c("Intercept", "Treaty created IGO", "Bureaucratic quality (lower)", "Water availability (lower)", "Treaty count", "Democratic dyad", "GDP p.c. (higher, logged)", "Power ratio", "Alliance", "Treaty created IGO $times$ Bureaucratic quality (lower)", "Treaty provides Monitoring", "Treaty provides Monitoring $times$ Bureaucratic quality (lower)", "Treaty provides Enforcement", "Treaty provides Enforcement $times$ Bureaucratic quality (lower)", "Treaty provides Conflict Resolution", "Treaty provides Conflict Resolution $times$ Bureaucratic quality (lower)"), 
               custom.note = "%stars. Standard errors in parentheses.", 
               digits = 3, 
               leading.zero = TRUE, 
               omit.coef = "factor", 
               reorder.coef = c(2, 3, 10, 11, 12, 13, 14, 15, 16, 4, 5, 6, 7, 8, 9, 1),
               center = TRUE, 
               caption = "Determinants of water-related cooperation.", 
               caption.above = TRUE, 
               label = "tab:robust2", 
               dcolumn = TRUE, 
               booktabs = TRUE, 
               use.packages = FALSE, 
               table = TRUE)
               
### Table 3: Seemingly unrelated regression  

# Note: be sure to correct N (n/2) and add correlation b/w the estimates by hand             
               
tab4 <- texreg(list(sur.bq),
		  file = "Output/tab_robust3.tex",
		  stars = 0.05,
		  # custom.model.names = c("Treaty institutionalization", "Cooperation"),
		  # custom.coef.names = c("Intercept", "Bureaucratic quality (lower)", "Water availability (lower)", "Treaty count", "Democratic dyad", "GDP p.c. (higher, logged)", "Power ratio", "Alliance", "Treaty institutionalization", "Treaty institutionalization $times$ Bureaucratic quality (lower)"),
		  custom.note = "%stars. ...",
		  digits = 3,
		  leading.zero = TRUE,
		  # reorder.coef = c(2, 9, 10, 3, 4, 5, 6, 7, 8, 1),
		  center = TRUE,
		  caption = "UPDATE N AND ADD CORRELATIONS!",
		  caption.above = TRUE,
		  label = "tab:sureg",
		  dcolumn = TRUE,
		  booktabs = TRUE,
		  use.packages = FALSE,
		  table = TRUE)               
               
               
#############################
### Plots
#############################

### 1.1 OLS BQ

icrg_burqual_min_std.sim <- sort(rep(quantile(ols.bq$x[, "icrg_burqual_min_std"], probs = c(0.1, 0.5, 0.9)), 5)) # length.out = 17 for more obs
instcoop_max2.sim <- rep(0:4, 3) # 17 accordingly
x.sim <- matrix(rep(as.numeric(apply(ols.bq$x[, c(4:9)], 2, function(x) quantile(x, probs = c(0.5)))), 5 * 3), ncol = 6, byrow = TRUE)

eff.ols.bq.dat <- as.matrix(cbind(1, instcoop_max2.sim, icrg_burqual_min_std.sim, x.sim, instcoop_max2.sim * icrg_burqual_min_std.sim))
eff.ols.bq.dat <- as.data.frame(eff.ols.bq.dat)
names(eff.ols.bq.dat) <- c(colnames(ols.bq$x))
rownames(eff.ols.bq.dat) <- NULL
pred <- predict(ols.bq, newdata = eff.ols.bq.dat, se.fit = TRUE, interval = "confidence", level = 0.90)

eff.ols.bq <- data.frame(eff.ols.bq.dat, pred$fit)
eff.ols.bq$icrg_burqual_min_std <- as.factor(eff.ols.bq$icrg_burqual_min_std)
levels(eff.ols.bq$icrg_burqual_min_std) <- c("BQ: 10th percentile", "BQ: Median", "BQ: 90th percentile")

# Three-in-one

p.ols.bq <- ggplot(dat = eff.ols.bq, aes(x = instcoop_max2, y = fit, color = icrg_burqual_min_std, group = icrg_burqual_min_std)) + geom_point(size = 5) + geom_ribbon(aes(x = instcoop_max2, ymin = lwr, ymax = upr, group = icrg_burqual_min_std, fill = icrg_burqual_min_std), alpha = 0.5, color = NA) + geom_line() + theme_jk() + labs(fill = "Bureaucratic \ncapacity") + theme(legend.position = c(0.25, 0.7)) + guides(color = FALSE, fill = FALSE) + xlab("Treaty institutionalization") + ylab("Predicted cooperation (continuous scale)") + ggtitle("") + annotate("text", x = 2.8, y = 0.7, label = "BQ: High", color = "#619CFF") + annotate("text", x = 2.5, y = 0.4, label = "BQ: Median", color = "#00BA38") + annotate("text", x = 1, y = 0.1, label = "BQ: Low", color = "#F8766D") # + scale_fill_brewer(palette="OrRd") + scale_color_brewer(palette="OrRd")

p.ols.bq_grayscale <- ggplot(dat = eff.ols.bq, aes(x = instcoop_max2, y = fit, color = icrg_burqual_min_std, group = icrg_burqual_min_std)) + geom_point(size = 5) + geom_ribbon(aes(x = instcoop_max2, ymin = lwr, ymax = upr, group = icrg_burqual_min_std, fill = icrg_burqual_min_std), alpha = 0.5, color = NA) + geom_line() + theme_jk() + labs(fill = "Bureaucratic \ncapacity") + theme(legend.position = c(0.25, 0.7)) + guides(color = FALSE, fill = FALSE) + xlab("Treaty institutionalization") + ylab("Predicted cooperation (continuous scale)") + ggtitle("") + scale_colour_grey(start = 0.8, end = 0) + scale_fill_grey(start = 0.8, end = 0) + annotate("text", x = 2.8, y = 0.7, label = "BQ: High", color = "gray0") + annotate("text", x = 2.5, y = 0.4, label = "BQ: Median", color = "gray30") + annotate("text", x = 1, y = 0.1, label = "BQ: Low", color = "gray60") 

ggsave(plot = p.ols.bq, filename = "Output/p_ols_bq.pdf", width = 4, height = 4)

# First differences (median minus low, etc.)

pred_for_fd <- predict(ols.bq, newdata = eff.ols.bq.dat, se.fit = TRUE)

pred_lo <- list()
pred_me <- list()
pred_hi <- list()

for(i in 1:5){
  pred_lo[[i]] <- rnorm(mean = pred_for_fd$fit[i], sd = pred_for_fd$se.fit[i], n = 100)
  pred_me[[i]] <- rnorm(mean = pred_for_fd$fit[i+5], sd = pred_for_fd$se.fit[i+5], n = 100)
  pred_hi[[i]] <- rnorm(mean = pred_for_fd$fit[i+10], sd = pred_for_fd$se.fit[i+10], n = 100)
}

pred_lo <- as.data.frame(pred_lo)
names(pred_lo) <- c("0", "1", "2", "3", "4")

pred_me <- as.data.frame(pred_me)
names(pred_me) <- c("0", "1", "2", "3", "4")

pred_hi <- as.data.frame(pred_hi)
names(pred_hi) <- c("0", "1", "2", "3", "4")

pred_fd_melo <- pred_me - pred_lo
pred_fd_hime <- pred_hi - pred_me

pred_fd_melo_sum <- data.frame(t(apply(X = pred_fd_melo, MARGIN = 2, FUN = function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))))
pred_fd_hime_sum <- data.frame(t(apply(X = pred_fd_hime, MARGIN = 2, FUN = function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))))

pred_fd_sum <- rbind(pred_fd_melo_sum, pred_fd_hime_sum)
pred_fd_sum$instcoop_max2 <- factor(rep(0:4, 2))
pred_fd_sum$Comparison <- factor(c(rep(1, 5), rep(2, 5)), labels = c("BQ: Median vs. Low", "BQ: High vs. Median"))

pred_fd_p <- ggplot(dat = pred_fd_sum, aes(x = instcoop_max2, y = X50.)) + 
  geom_segment(aes(xend = instcoop_max2, y = X5., yend = X95.)) + 
  geom_point() + theme_jk()  +
  facet_wrap(~ Comparison, ncol = 1) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Treaty institutionalization") + 
  ylab("Difference in\npredicted cooperation (continuous scale)") + 
  ggtitle("")

ggsave(pred_fd_p, file = "Output/fd_ols_bq.pdf", width = 4, height = 4)

pdf("Output/pfd_ols_bq.pdf", width = 4, height = 6)
grid.arrange(p.ols.bq, pred_fd_p, ncol = 1, heights = c(1.5, 1))
dev.off()


# Marginal effects

ols.bq.mep <- ggintfun(ols.bq, 
                       varnames = c("instcoop_max2", "icrg_burqual_min_std"), 
                       varlabs = c("Treaty institutionalization", "Bureaucratic quality"), 
                       rug = TRUE, jitter_factor = 5) + 
  theme_jk() + 
  scale_y_continuous(limits = c(-0.05, 0.25))
ggsave(plot = ols.bq.mep, filename = "Output/megg_ols_bq.pdf", width = 4, height = 4)

ols.inst.mep <- ggintfun(ols.bq, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = TRUE, jitter_factor = 1) + theme_jk()
ggsave(plot = ols.inst.mep, filename = "Output/megg_ols_inst.pdf", width = 4, height = 4)

# Figure 1

pdf("Output/fig1_color.pdf", width = 11, height = 4.5)
grid.arrange(ols.bq.mep + ggtitle("Conditional effect"), p.ols.bq + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 3, widths = c(1, 1, 0.6))
dev.off()

pdf("Output/fig1.pdf", width = 11, height = 4.5)
grid.arrange(ols.bq.mep + ggtitle("Conditional effect"), p.ols.bq_grayscale + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 3, widths = c(1, 1, 0.6))
dev.off()

tiff("Output/fig1.tif", res = 350, compression = "lzw", width = 11, height = 4.5, units = "in")
grid.arrange(ols.bq.mep + ggtitle("Conditional effect"), p.ols.bq_grayscale + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 3, widths = c(1, 1, 0.6))
dev.off()

###  1.2 Probit BQ

# Three-in-one

# Simulated data

icrg_burqual_min_std.sim <- sort(rep(quantile(ols.bq$x[, "icrg_burqual_min_std"], probs = c(0.1, 0.5, 0.9)), 5)) # length.out = 17 for more obs
instcoop_max2.sim <- rep(0:4, 3) # 17 accordingly
x.sim <- matrix(rep(as.numeric(apply(ols.bq$x[, c(4:9)], 2, function(x) quantile(x, probs = c(0.5)))), 5 * 3), ncol = 6, byrow = TRUE)

eff.probit.bq <- effect(term = "instcoop_max2 * icrg_burqual_min_std", mod = probit.bq, xlevels = list(instcoop_max2 = c(0, 1, 2, 3, 4), icrg_burqual_min_std = quantile(probit.bq$x[, "icrg_burqual_min_std"], probs = c(0.1, 0.5, 0.9))), se = TRUE, confidence.level = 0.90, typical = median)

eff.probit.bq.dat <- cbind(eff.probit.bq$x, pnorm(eff.probit.bq$fit), pnorm(eff.probit.bq$lower), pnorm(eff.probit.bq$upper))
names(eff.probit.bq.dat) <- c("instcoop_max2", "icrg_burqual_min_std", "fit", "lwr", "upr")
eff.probit.bq.dat$icrg_burqual_min_std <- as.factor(eff.probit.bq.dat$icrg_burqual_min_std)
levels(eff.probit.bq.dat$icrg_burqual_min_std) <- c("BQ: 10th percentile", "BQ: Median", "BQ: 90th percentile")

p.probit.bq <- ggplot(dat = eff.probit.bq.dat, aes(x = instcoop_max2, y = fit, color = icrg_burqual_min_std, group = icrg_burqual_min_std)) + geom_point(size = 5) + geom_ribbon(aes(x = instcoop_max2, ymin = lwr, ymax = upr, group = icrg_burqual_min_std, fill = icrg_burqual_min_std), alpha = 0.5, color = NA) + geom_line() + theme_jk() + labs(fill = "Bureaucratic \ncapacity") + theme(legend.position = c(0.25, 0.7)) + guides(color = FALSE, fill = FALSE) + xlab("Treaty institutionalization") + ylab("Predicted probability of 1+ cooperative event") + ggtitle("") + ylim(c(0, 0.33)) + annotate("text", x = 2, y = 0.2, label = "BQ: High", color = "#619CFF") + annotate("text", x = 1.5, y = 0.075, label = "BQ: Median", color = "#00BA38") + annotate("text", x = 0.8, y = 0.05, label = "BQ: Low", color = "#F8766D")

p.probit.bq_grayscale <- ggplot(dat = eff.probit.bq.dat, aes(x = instcoop_max2, y = fit, color = icrg_burqual_min_std, group = icrg_burqual_min_std)) + geom_point(size = 5) + geom_ribbon(aes(x = instcoop_max2, ymin = lwr, ymax = upr, group = icrg_burqual_min_std, fill = icrg_burqual_min_std), alpha = 0.5, color = NA) + geom_line() + theme_jk() + labs(fill = "Bureaucratic \ncapacity") + theme(legend.position = c(0.25, 0.7)) + guides(color = FALSE, fill = FALSE) + xlab("Treaty institutionalization") + ylab("Predicted probability of 1+ cooperative event") + ggtitle("") + ylim(c(0, 0.33)) + scale_colour_grey(start = 0.8, end = 0) + scale_fill_grey(start = 0.8, end = 0) + annotate("text", x = 2, y = 0.2, label = "BQ: High", color = "gray0") + annotate("text", x = 1.5, y = 0.075, label = "BQ: Median", color = "gray30") + annotate("text", x = 0.8, y = 0.05, label = "BQ: Low", color = "gray60")

ggsave(plot = p.probit.bq, filename = "Output/p_probit_bq.pdf", width = 4, height = 4)

# First differences (median minus low, etc.)

pred_for_fd <- predict(probit.bq, newdata = eff.ols.bq.dat, type = "response", se.fit = TRUE)

pred_lo <- list()
pred_me <- list()
pred_hi <- list()

for(i in 1:5){
  pred_lo[[i]] <- rnorm(mean = pred_for_fd$fit[i], sd = pred_for_fd$se.fit[i], n = 100)
  pred_me[[i]] <- rnorm(mean = pred_for_fd$fit[i+5], sd = pred_for_fd$se.fit[i+5], n = 100)
  pred_hi[[i]] <- rnorm(mean = pred_for_fd$fit[i+10], sd = pred_for_fd$se.fit[i+10], n = 100)
}

pred_lo <- as.data.frame(pred_lo)
names(pred_lo) <- c("0", "1", "2", "3", "4")

pred_me <- as.data.frame(pred_me)
names(pred_me) <- c("0", "1", "2", "3", "4")

pred_hi <- as.data.frame(pred_hi)
names(pred_hi) <- c("0", "1", "2", "3", "4")

pred_fd_melo <- pred_me - pred_lo
pred_fd_hime <- pred_hi - pred_me

pred_fd_melo_sum <- data.frame(t(apply(X = pred_fd_melo, MARGIN = 2, FUN = function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))))
pred_fd_hime_sum <- data.frame(t(apply(X = pred_fd_hime, MARGIN = 2, FUN = function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))))

pred_fd_sum <- rbind(pred_fd_melo_sum, pred_fd_hime_sum)
pred_fd_sum$instcoop_max2 <- factor(rep(0:4, 2))
pred_fd_sum$Comparison <- factor(c(rep(1, 5), rep(2, 5)), labels = c("BQ: Median vs. Low", "BQ: High vs. Median"))

pred_fd_p <- ggplot(dat = pred_fd_sum, aes(x = instcoop_max2, y = X50.)) + 
  geom_segment(aes(xend = instcoop_max2, y = X5., yend = X95.)) + 
  geom_point() + theme_jk()  +
  facet_wrap(~ Comparison, ncol = 1) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Treaty institutionalization") + 
  ylab("Difference in\nPr(1+ cooperative event)") + 
  ggtitle("")

ggsave(pred_fd_p, file = "Output/fd_probit_bq.pdf", width = 4, height = 4)

pdf("Output/pfd_probit_bq.pdf", width = 4, height = 6)
grid.arrange(p.probit.bq, pred_fd_p, ncol = 1, heights = c(1.5, 1))
dev.off()

# Marginal effects
# CIs Via sampling from MVN 

# simulate coefficients based on the model

beta <- coef(probit.bq)
covvar <- summary(probit.bq)$cov.unscaled
coefs.sim <- MASS::mvrnorm(n = 10000, mu = beta, Sigma = covvar)

# Simulate data for different values of moderator
# Moderator: bureaucratic quality

sim_Intercept <- 1
sim_instcoop_max2 <- 0
sim_icrg_burqual_min_std <- seq(0, 1, by = 0.1)
sim_wres_min_log <- median(probit.bq$model$wres_min_log)
sim_ntreaties2 <- median(probit.bq$model$ntreaties2)
sim_dyaddem7 <- median(probit.bq$model$dyaddem7)
sim_gdpmaxPWT_log <- median(probit.bq$model$gdpmaxPWT_log)
sim_lnpowerratio <- median(probit.bq$model$lnpowerratio)
sim_alliance <- median(probit.bq$model$alliance)
sim_instcoop_max2_icrg_burqual_min_std <- sim_instcoop_max2 * sim_icrg_burqual_min_std 

fake_inst0 <- cbind(sim_Intercept, sim_instcoop_max2, sim_icrg_burqual_min_std, sim_wres_min_log, sim_ntreaties2, sim_dyaddem7, sim_gdpmaxPWT_log, sim_lnpowerratio, sim_alliance, sim_instcoop_max2_icrg_burqual_min_std)

sim_instcoop_max2 <- 1
sim_instcoop_max2_icrg_burqual_min_std <- sim_instcoop_max2 * sim_icrg_burqual_min_std
fake_inst1 <- cbind(sim_Intercept, sim_instcoop_max2, sim_icrg_burqual_min_std, sim_wres_min_log, sim_ntreaties2, sim_dyaddem7, sim_gdpmaxPWT_log, sim_lnpowerratio, sim_alliance, sim_instcoop_max2_icrg_burqual_min_std)

sim_instcoop_max2 <- 2
sim_instcoop_max2_icrg_burqual_min_std <- sim_instcoop_max2 * sim_icrg_burqual_min_std
fake_inst2 <- cbind(sim_Intercept, sim_instcoop_max2, sim_icrg_burqual_min_std, sim_wres_min_log, sim_ntreaties2, sim_dyaddem7, sim_gdpmaxPWT_log, sim_lnpowerratio, sim_alliance, sim_instcoop_max2_icrg_burqual_min_std)

sim_instcoop_max2 <- 3
sim_instcoop_max2_icrg_burqual_min_std <- sim_instcoop_max2 * sim_icrg_burqual_min_std
fake_inst3 <- cbind(sim_Intercept, sim_instcoop_max2, sim_icrg_burqual_min_std, sim_wres_min_log, sim_ntreaties2, sim_dyaddem7, sim_gdpmaxPWT_log, sim_lnpowerratio, sim_alliance, sim_instcoop_max2_icrg_burqual_min_std)

sim_instcoop_max2 <- 4
sim_instcoop_max2_icrg_burqual_min_std <- sim_instcoop_max2 * sim_icrg_burqual_min_std
fake_inst4 <- cbind(sim_Intercept, sim_instcoop_max2, sim_icrg_burqual_min_std, sim_wres_min_log, sim_ntreaties2, sim_dyaddem7, sim_gdpmaxPWT_log, sim_lnpowerratio, sim_alliance, sim_instcoop_max2_icrg_burqual_min_std)

# multiply the fake data matrix with the matrix of simulated coefficients to obtain marginal effects
results0 <- coefs.sim %*% t(fake_inst0)
results1 <- coefs.sim %*% t(fake_inst1)
results2 <- coefs.sim %*% t(fake_inst2)
results3 <- coefs.sim %*% t(fake_inst3)
results4 <- coefs.sim %*% t(fake_inst4)

# evaluate each marginal effect on the standardized normal CDF, to obtain predicted probabilities
phat0 <- apply(results0, 1:2, function(x) pnorm(x))
phat1 <- apply(results1, 1:2, function(x) pnorm(x))
phat2 <- apply(results2, 1:2, function(x) pnorm(x))
phat3 <- apply(results3, 1:2, function(x) pnorm(x))
phat4 <- apply(results4, 1:2, function(x) pnorm(x))

# subtract the estimated probabilities from each other
fd01 <- phat1 - phat0
fd12 <- phat2 - phat1
fd23 <- phat3 - phat2
fd34 <- phat4 - phat3

# What's the relevant change?
quantile(probit.bq$model$instcoop_max2, probs = c(0.25, 0.75))
# 0 to 2
fd02 <- phat2 - phat0
quantile(probit.bq$model$instcoop_max2, probs = c(0.1, 0.9))
# 0 to 3
fd03 <- phat3 - phat0
# 0 to 4
fd04 <- phat4 - phat0

# graph the results
fd01_means <- apply(fd01, 2, mean)
fd01_sdUpper <- apply(fd01, 2, function(x) quantile(x, .95))
fd01_sdLower <- apply(fd01, 2, function(x) quantile(x, .05))
fd02_means <- apply(fd02, 2, mean)
fd02_sdUpper <- apply(fd02, 2, function(x) quantile(x, .95))
fd02_sdLower <- apply(fd02, 2, function(x) quantile(x, .05))
fd03_means <- apply(fd03, 2, mean)
fd03_sdUpper <- apply(fd03, 2, function(x) quantile(x, .95))
fd03_sdLower <- apply(fd03, 2, function(x) quantile(x, .05))
fd04_means <- apply(fd04, 2, mean)
fd04_sdUpper <- apply(fd04, 2, function(x) quantile(x, .95))
fd04_sdLower <- apply(fd04, 2, function(x) quantile(x, .05))

plot(x = sim_icrg_burqual_min_std, y = fd03_means, type = "l", ylim = c(0, 0.2))
lines(x = sim_icrg_burqual_min_std, y = fd03_sdUpper, lty = 2)
lines(x = sim_icrg_burqual_min_std, y = fd03_sdLower, lty = 2)            
               
probit.bq_fd <- data.frame(fd.mean = c(fd01_means, fd02_means, fd03_means, fd04_means), 	
					 fd.lo = c(fd01_sdLower, fd02_sdLower, fd03_sdLower, fd04_sdLower),
					 fd.hi = c(fd01_sdUpper, fd02_sdUpper, fd03_sdUpper, fd04_sdUpper),
					 icrg_burqual_min_std = rep(sim_icrg_burqual_min_std, 4),
					 instcoop_change = c(rep("0 -> 1", 11), rep("0 -> 2", 11), rep("0 -> 3", 11), rep("0 -> 4", 11)))

probit.bq.me3p <- ggplot(data = probit.bq_fd[probit.bq_fd$instcoop_change == "0 -> 3", ], aes(x = icrg_burqual_min_std, y = fd.mean)) +
          geom_hline(yintercept = 0, linetype = "dashed") +
          geom_ribbon(aes(x = icrg_burqual_min_std, ymin = fd.lo, ymax = fd.hi), alpha = 0.25, color = NA) +
          geom_line() + 
          # geom_line(aes(x = icrg_burqual_min_std, y = fd.lo), linetype = "dashed") + 
          # geom_line(aes(x = icrg_burqual_min_std, y = fd.hi), linetype = "dashed") +
          xlab("Bureaucratic quality") + ylab("Effect of treaty institutionalization\non Pr(1+ cooperative event)") +
          theme_jk() + scale_y_continuous(limits = c(-0.05, 0.25))

rug.dat <- data.frame(probit.bq$model["icrg_burqual_min_std"])
probit.bq.me3p <- probit.bq.me3p + geom_rug(data = rug.dat, aes(x = jitter(rug.dat[, 1], factor = 5), y = 0), sides = "b", size = 0.1)

ggsave(plot = probit.bq.me3p, filename = "Output/megg_probit_bq.pdf", width = 4, height = 4)

# Figure 2

pdf("Output/fig2_color.pdf", width = 11, height = 4.5)
grid.arrange(probit.bq.me3p + ggtitle("Conditional effect"), p.probit.bq + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 3, widths = c(1, 1, 0.6))
dev.off()

pdf("Output/fig2.pdf", width = 11, height = 4.5)
grid.arrange(probit.bq.me3p + ggtitle("Conditional effect"), p.probit.bq_grayscale + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 3, widths = c(1, 1, 0.6))
dev.off()

tiff("Output/fig2.tif", res = 350, compression = "lzw", width = 11, height = 4.5, units = "in")
grid.arrange(probit.bq.me3p + ggtitle("Conditional effect"), p.probit.bq_grayscale + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 3, widths = c(1, 1, 0.6))
dev.off()

# Simulate data for different values of moderator
# Moderator: treaty institutionalization

sim_Intercept <- 1
sim_instcoop_max2 <- c(1, 2, 3, 4)
sim_icrg_burqual_min_std <- 0
sim_wres_min_log <- median(probit.bq$model$wres_min_log)
sim_ntreaties2 <- median(probit.bq$model$ntreaties2)
sim_dyaddem7 <- median(probit.bq$model$dyaddem7)
sim_gdpmaxPWT_log <- median(probit.bq$model$gdpmaxPWT_log)
sim_lnpowerratio <- median(probit.bq$model$lnpowerratio)
sim_alliance <- median(probit.bq$model$alliance)
sim_instcoop_max2_icrg_burqual_min_std <- sim_instcoop_max2 * sim_icrg_burqual_min_std 

fake_bqlo <- cbind(sim_Intercept, sim_instcoop_max2, sim_icrg_burqual_min_std, sim_wres_min_log, sim_ntreaties2, sim_dyaddem7, sim_gdpmaxPWT_log, sim_lnpowerratio, sim_alliance, sim_instcoop_max2_icrg_burqual_min_std)

sim_icrg_burqual_min_std <- 0.83
sim_instcoop_max2_icrg_burqual_min_std <- sim_instcoop_max2 * sim_icrg_burqual_min_std
fake_bqhi <- cbind(sim_Intercept, sim_instcoop_max2, sim_icrg_burqual_min_std, sim_wres_min_log, sim_ntreaties2, sim_dyaddem7, sim_gdpmaxPWT_log, sim_lnpowerratio, sim_alliance, sim_instcoop_max2_icrg_burqual_min_std)

# multiply the fake data matrix with the matrix of simulated coefficients to obtain marginal effects
resultslo <- coefs.sim %*% t(fake_bqlo)
resultshi <- coefs.sim %*% t(fake_bqhi)

# evaluate each marginal effect on the standardized normal CDF, to obtain predicted probabilities
phatlo <- apply(resultslo, 1:2, function(x) pnorm(x))
phathi <- apply(resultshi, 1:2, function(x) pnorm(x))

# subtract the estimated probabilities from each other
fdhilo <- phathi - phatlo

# What's the relevant change?
quantile(probit.bq$model$icrg_burqual_min_std, probs = c(0.1, 0.9))
# 0 to 0.83

# graph the results
fdhilo_means <- apply(fdhilo, 2, mean)
fdhilo_sdUpper <- apply(fdhilo, 2, function(x) quantile(x, .95))
fdhilo_sdLower <- apply(fdhilo, 2, function(x) quantile(x, .05))

probit.inst_fd <- data.frame(fd.mean = fdhilo_means, 	
                             fd.lo = fdhilo_sdLower,
                             fd.hi = fdhilo_sdUpper,
                             instcoop_max2 = sim_instcoop_max2,
                             instcoop_change = c(rep("Low -> High", 4)))

probit.inst.me3p <- ggplot(data = probit.inst_fd, aes(x = instcoop_max2, y = fd.mean)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_ribbon(aes(x = instcoop_max2, ymin = fd.lo, ymax = fd.hi), alpha = 0.25, color = NA) +
  geom_line() + 
  xlab("Treaty institutionalization") + ylab("Effect of bureaucratic quality\non Pr(Cooperation)") +
  theme_jk() 

ggsave(plot = probit.inst.me3p, filename = "Output/megg_probit_inst.pdf", width = 4, height = 4)

###  2. Years with events only

# Marginal effects

ols.bq.eo.mep <- ggintfun(ols.bq.eo, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = TRUE, twoways = FALSE, jitter_factor = 5)
ggsave(plot = ols.bq.eo.mep, filename = "Output/megg_ols_bq_eo.pdf", width = 4, height = 4)

ols.inst.eo.mep <- ggintfun(ols.bq.eo, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = TRUE, twoways = FALSE, jitter_factor = 1)
ggsave(plot = ols.inst.eo.mep, filename = "Output/megg_ols_inst_eo.pdf", width = 4, height = 4)


###  3. Median cooperation instead of mean/average

# Marginal effects

ols.bq.eo.mep <- ggintfun(ols.bq.median, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE, twoways = FALSE, jitter_factor = 5) + scale_y_continuous(limits = c(-0.05, 0.25))
ggsave(plot = ols.bq.eo.mep, filename = "Output/megg_ols_bq_median.pdf", width = 4, height = 4)

ols.inst.eo.mep <- ggintfun(ols.bq.median, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE, twoways = FALSE, jitter_factor = 1)
ggsave(plot = ols.inst.eo.mep, filename = "Output/megg_ols_inst_median.pdf", width = 4, height = 4)

###  4. Indicator for no treaty
# ols.bq.3w

icrg_burqual_min_std.sim <- sort(rep(quantile(ols.bq.3w$x[, "icrg_burqual_min_std"], probs = c(0.1, 0.5, 0.9)), 6)) # length.out = 18 for more obs
instcoop_max2.sim <- c(0, 0, 1, 2, 3, 4, 0, 0, 1, 2, 3, 4, 0, 0, 1, 2, 3, 4) # 18 accordingly
notreaties.sim <- c(1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)

x.sim <- matrix(rep(as.numeric(apply(ols.bq.3w$x[, c(5:10)], 2, function(x) quantile(x, probs = c(0.5)))), 6 * 3), ncol = 6, byrow = TRUE)

eff.ols.bq.3w.dat <- as.matrix(cbind(1, instcoop_max2.sim, icrg_burqual_min_std.sim, notreaties.sim, x.sim, instcoop_max2.sim * icrg_burqual_min_std.sim, instcoop_max2.sim * notreaties.sim, icrg_burqual_min_std.sim * notreaties.sim, instcoop_max2.sim * icrg_burqual_min_std.sim * notreaties.sim))
eff.ols.bq.3w.dat <- as.data.frame(eff.ols.bq.3w.dat)
names(eff.ols.bq.3w.dat) <- c(colnames(ols.bq.3w$x))
rownames(eff.ols.bq.3w.dat) <- NULL
pred <- predict(ols.bq.3w, newdata = eff.ols.bq.3w.dat, se.fit = TRUE, interval = "confidence", level = 0.90)

eff.ols.bq.3w <- data.frame(eff.ols.bq.3w.dat, pred$fit)
eff.ols.bq.3w$icrg_burqual_min_std <- as.factor(eff.ols.bq.3w$icrg_burqual_min_std)
levels(eff.ols.bq.3w$icrg_burqual_min_std) <- c("BQ: 10th percentile", "BQ: Median", "BQ: 90th percentile")

# Three-in-one

# First, create two datasets, one for notreaties == 1, and one for notreaties == 0

eff.ols.bq.3w.notreaties <- eff.ols.bq.3w[eff.ols.bq.3w$notreaties == 1, ]
eff.ols.bq.3w.notreaties$instcoop_max2 <- c(-0.15, 0, 0.15)
eff.ols.bq.3w.treaties <- eff.ols.bq.3w[eff.ols.bq.3w$notreaties == 0, ]

p.ols.bq.3w <- ggplot(dat = eff.ols.bq.3w.treaties, aes(x = instcoop_max2, y = fit, color = icrg_burqual_min_std, group = icrg_burqual_min_std)) + geom_point(size = 5) +
  geom_ribbon(aes(x = instcoop_max2, ymin = lwr, ymax = upr, group = icrg_burqual_min_std, fill = icrg_burqual_min_std), alpha = 0.5, color = NA) + 
  geom_line() + 
  geom_segment(dat = eff.ols.bq.3w.notreaties, aes(x = instcoop_max2 - 1, xend = instcoop_max2 - 1, y = lwr, yend = upr, color = icrg_burqual_min_std), size = 1.5, alpha = 0.5) +
  geom_point(dat = eff.ols.bq.3w.notreaties, aes(x = instcoop_max2 - 1, y = fit, color = icrg_burqual_min_std), size = 5) +
  theme_jk()  + 
  scale_x_continuous(breaks = c(-1:4), labels = c("No treaty", "0", "1", "2", "3", "4")) +
  labs(fill = "Bureaucratic \ncapacity") + theme(legend.position = c(0.25, 0.7)) + 
  guides(color = FALSE, fill = FALSE) + 
  xlab("Treaty institutionalization") + 
  ylab("Predicted cooperation (continuous scale)") + 
  annotate("text", x = 1.25, y = 0.6, label = "BQ: High", colour = "#619CFF") + 
  annotate("text", x = 2.5, y = 0.35, label = "BQ: Median", colour = "#00BA38") + 
  annotate("text", x = 1, y = 0, label = "BQ: Low", colour = "#F8766D") + 
  annotate("text", x = -1, y = 0.075, label = "BQ: High", colour = "#619CFF", hjust = 0.2) + 
  annotate("text", x = -0.85, y = 0.3, label = "BQ:\nMedian", colour = "#00BA38", hjust = 0) + 
  annotate("text", x = -1.2, y = 0.4, label = "BQ: Low", colour = "#F8766D", hjust = 0) 

p.ols.bq.3w_grayscale <- ggplot(dat = eff.ols.bq.3w.treaties, aes(x = instcoop_max2, y = fit, color = icrg_burqual_min_std, group = icrg_burqual_min_std)) + geom_point(size = 5) +
  geom_ribbon(aes(x = instcoop_max2, ymin = lwr, ymax = upr, group = icrg_burqual_min_std, fill = icrg_burqual_min_std), alpha = 0.5, color = NA) + 
  geom_line() + 
  geom_segment(dat = eff.ols.bq.3w.notreaties, aes(x = instcoop_max2 - 1, xend = instcoop_max2 - 1, y = lwr, yend = upr, color = icrg_burqual_min_std), size = 1.5, alpha = 0.5) +
  geom_point(dat = eff.ols.bq.3w.notreaties, aes(x = instcoop_max2 - 1, y = fit, color = icrg_burqual_min_std), size = 5) +
  theme_jk()  + 
  scale_x_continuous(breaks = c(-1:4), labels = c("No treaty", "0", "1", "2", "3", "4")) +
  labs(fill = "Bureaucratic \ncapacity") + theme(legend.position = c(0.25, 0.7)) + 
  guides(color = FALSE, fill = FALSE) + 
  xlab("Treaty institutionalization") + 
  ylab("Predicted cooperation (continuous scale)") + 
  annotate("text", x = 1.25, y = 0.6, label = "BQ: High", color = "gray0") + 
  annotate("text", x = 2.5, y = 0.35, label = "BQ: Median", color = "gray30") + 
  annotate("text", x = 1, y = 0, label = "BQ: Low", color = "gray60") + 
  annotate("text", x = -1, y = 0.075, label = "BQ: High", color = "gray0", hjust = 0.2) + 
  annotate("text", x = -0.85, y = 0.3, label = "BQ:\nMedian", color = "gray30", hjust = 0) + 
  annotate("text", x = -1.2, y = 0.4, label = "BQ: Low", color = "gray60", hjust = 0) + 
  scale_colour_grey(start = 0.8, end = 0) + scale_fill_grey(start = 0.8, end = 0)

ggsave(plot = p.ols.bq.3w, filename = "Output/p_ols_bq_3w.pdf", width = 5, height = 4)

# First differences (median minus low, etc.)

pred_for_fd <- predict(ols.bq.3w, newdata = eff.ols.bq.3w.dat, se.fit = TRUE)

pred_lo <- list()
pred_me <- list()
pred_hi <- list()

for(i in 1:6){
  pred_lo[[i]] <- rnorm(mean = pred_for_fd$fit[i], sd = pred_for_fd$se.fit[i], n = 100)
  pred_me[[i]] <- rnorm(mean = pred_for_fd$fit[i+6], sd = pred_for_fd$se.fit[i+5], n = 100)
  pred_hi[[i]] <- rnorm(mean = pred_for_fd$fit[i+12], sd = pred_for_fd$se.fit[i+10], n = 100)
}

pred_lo <- as.data.frame(pred_lo)
names(pred_lo) <- c("No treaties", "0", "1", "2", "3", "4")

pred_me <- as.data.frame(pred_me)
names(pred_me) <- c("No treaties", "0", "1", "2", "3", "4")

pred_hi <- as.data.frame(pred_hi)
names(pred_hi) <- c("No treaties", "0", "1", "2", "3", "4")

pred_fd_melo <- pred_me - pred_lo
pred_fd_hime <- pred_hi - pred_me

pred_fd_melo_sum <- data.frame(t(apply(X = pred_fd_melo, MARGIN = 2, FUN = function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))))
pred_fd_hime_sum <- data.frame(t(apply(X = pred_fd_hime, MARGIN = 2, FUN = function(x) quantile(x, probs = c(0.05, 0.5, 0.95)))))

pred_fd_sum <- rbind(pred_fd_melo_sum, pred_fd_hime_sum)
pred_fd_sum$instcoop_max2 <- factor(rep(0:5, 2), labels =  c("No\ntreaty", "0", "1", "2", "3", "4"))
pred_fd_sum$Comparison <- factor(c(rep(1, 6), rep(2, 6)), labels = c("BQ: Median vs. Low", "BQ: High vs. Median"))

pred_fd_p <- ggplot(dat = pred_fd_sum, aes(x = instcoop_max2, y = X50.)) + 
  geom_segment(aes(xend = instcoop_max2, y = X5., yend = X95.)) + 
  geom_point() + theme_jk()  +
  facet_wrap(~ Comparison, ncol = 1) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  xlab("Treaty institutionalization") + 
  ylab("Difference in predicted cooperation\n(continuous scale)") + 
  theme(axis.text.x = element_text(size = 8)) +
  ggtitle("")

ggsave(pred_fd_p, file = "Output/fd_ols_bq_3w.pdf", width = 4, height = 4)

pdf("Output/pfd_ols_bq_3w.pdf", width = 5, height = 6)
grid.arrange(p.ols.bq.3w, pred_fd_p, ncol = 1, heights = c(1.5, 1))
dev.off()

# Figure 3

pdf("Output/fig3_color.pdf", width = 7, height = 4.5)
grid.arrange(p.ols.bq.3w + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 2, widths = c(1, 0.6))
dev.off()

pdf("Output/fig3.pdf", width = 7, height = 4.5)
grid.arrange(p.ols.bq.3w_grayscale + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 2, widths = c(1, 0.6))
dev.off()

tiff("Output/fig3.tif", res = 350, compression = "lzw", width = 7, height = 4.5, units = "in")
grid.arrange(p.ols.bq.3w_grayscale + ggtitle("Predicted cooperation"), pred_fd_p + ggtitle("Differences"), ncol = 2, widths = c(1, 0.6))
dev.off()

# Marginal effects

icrg_burqual_min_std.range <- c(seq(0, 1, by = 0.01))
eff <- coef(ols.bq.3w)[2] + coef(ols.bq.3w)[11] * icrg_burqual_min_std.range 
eff.var <- vcov(ols.bq.3w)[2, 2] + (icrg_burqual_min_std.range^2 * vcov(ols.bq.3w)[11, 11]) + (2 * icrg_burqual_min_std.range * vcov(ols.bq.3w)[2, 11])
eff.se <- sqrt(eff.var)
low <- eff - qt(0.975, ols.bq.3w$df.residual) * eff.se
up <- eff + qt(0.975, ols.bq.3w$df.residual) * eff.se

rug.dat <- data.frame(ols.bq.3w$model["icrg_burqual_min_std"])

plot.dat <- data.frame(icrg_burqual_min_std.range, eff, eff.se, low, up)

megg.ols.bq.3w <- ggplot(data = plot.dat, aes(x = icrg_burqual_min_std.range, y = eff), environment = environment()) +     # env. important for rug
  geom_hline(yintercept = 0, linetype = "dashed") + theme_jk() +
  geom_ribbon(aes(x = icrg_burqual_min_std.range, ymin = low, ymax = up), alpha = 0.25, color = NA) +
  geom_line() + 
  scale_y_continuous(limits = c(-0.07, 0.25)) + #, breaks = c(-0.05, 0, 0.05, 0.1, 0.15, 0.2, 0.25)) +
  xlab("Bureaucratic quality") + ylab("Effect of Treaty institutionalization") # +

ggsave(plot = megg.ols.bq.3w, filename = "Output/megg_ols_bq_3w.pdf", width = 4, height = 4)


inst_coop_max2.sim <- c(0, 1, 2, 3, 4)
eff <- coef(ols.bq.3w)[3] + coef(ols.bq.3w)[11] * inst_coop_max2.sim
eff.var <- vcov(ols.bq.3w)[3, 3] + (inst_coop_max2.sim^2 * vcov(ols.bq.3w)[11, 11]) + (2 * inst_coop_max2.sim * vcov(ols.bq.3w)[3, 11])
eff.se <- sqrt(eff.var)
low <- eff - qt(0.975, ols.bq.3w$df.residual) * eff.se
up <- eff + qt(0.975, ols.bq.3w$df.residual) * eff.se

rug.dat <- data.frame(ols.bq.3w$model["instcoop_max2"])

plot.dat <- data.frame(treaty = c(1:5), eff, eff.se, low, up)

eff <- coef(ols.bq.3w)[3] + coef(ols.bq.3w)[13] * 1
eff.var <- vcov(ols.bq.3w)[3, 3] + (1 * vcov(ols.bq.3w)[12, 12]) + (2 * 1 * vcov(ols.bq.3w)[3, 12])
eff.se <- sqrt(eff.var)
low <- eff - qt(0.975, ols.bq.3w$df.residual) * eff.se
up <- eff + qt(0.975, ols.bq.3w$df.residual) * eff.se

plot.dat <- rbind(plot.dat, c(0, eff, eff.se, low, up))
plot.dat <- arrange(plot.dat, treaty)

plot.dat$treaty <- factor(plot.dat$treaty, labels = c("No\ntreaties", "0", "1", "2", "3", "4"))

megg.ols.inst.3w <- ggplot(data = plot.dat, aes(x = treaty, y = eff)) + 
  geom_hline(yintercept = 0, linetype = "dashed") + theme_jk() +
  geom_segment(aes(x = treaty, xend = treaty, y = low, yend = up)) +
  geom_point() +
  xlab("Treaty institutionalization") + ylab("Effect of Bureaucratic Quality") # +

ggsave(plot = megg.ols.inst.3w, filename = "Output/megg_ols_inst_3w.pdf", width = 4, height = 4)


###  5. Dyad FE
# ols.bq.dfe

ols.bq.dfe.mep <- ggintfun(ols.bq.dfe, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE) + scale_y_continuous(limits = c(NA, 0.25))
ggsave(plot = ols.bq.dfe.mep, filename = "Output/megg_ols_bq_dfe.pdf", width = 4, height = 4)

ols.inst.dfe.mep <- ggintfun(ols.bq.dfe, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE)
ggsave(plot = ols.inst.dfe.mep, filename = "Output/megg_ols_inst_dfe.pdf", width = 4, height = 4)


###  6a. Dyad FELDV
# ols.bq.dfeldv

ols.bq.dfeldv.mep <- ggintfun(ols.bq.dfeldv, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE) + scale_y_continuous(limits = c(NA, 0.25))
ggsave(plot = ols.bq.dfeldv.mep, filename = "Output/megg_ols_bq_dfeldv.pdf", width = 4, height = 4)

ols.inst.dfeldv.mep <- ggintfun(ols.bq.dfeldv, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE)
ggsave(plot = ols.inst.dfeldv.mep, filename = "Output/megg_ols_inst_dfeldv.pdf", width = 4, height = 4)

###  6b. Arellano-Bond
# ols.abgmm

# ols.abgmm$coefficients
# vcovHC(ols.abgmm)

icrg_burqual_min_std <- seq(0, 1, by = 0.05)
b1.instcoop <- ols.abgmm$coefficients[1]
b3.interact <- ols.abgmm$coefficients[9]
var1.instcoop <- plm::vcovHC(ols.abgmm)[1, 1]
var3.interact <- plm::vcovHC(ols.abgmm)[9, 9]
cov13.instcoop_interact <- plm::vcovHC(ols.abgmm)[1, 9]
mf.instcoop <- b1.instcoop + b3.interact * icrg_burqual_min_std
var.instcoop <- var1.instcoop + icrg_burqual_min_std^2 * var3.interact + 2 * icrg_burqual_min_std * cov13.instcoop_interact
se.instcoop <- sqrt(var.instcoop)
low1.instcoop <- mf.instcoop - 1.96 * se.instcoop
up1.instcoop <- mf.instcoop + 1.96 * se.instcoop

ols.abgmm.plot.dat <- data.frame(s2 = icrg_burqual_min_std, eff1 = mf.instcoop, se.eff1 = se.instcoop, low1 = low1.instcoop, up1 = up1.instcoop)

ols.bq.abgmm.mep <- ggplot(data = ols.abgmm.plot.dat, aes(x = s2, y = eff1)) + 
  geom_hline(yintercept = 0, linetype = "dashed") + theme_jk() +
  geom_ribbon(aes(x = s2, ymin = low1, ymax = up1), alpha = 0.25, color = NA) +
  geom_line() + 
  scale_x_continuous(breaks = c(seq(0, 1, by = 0.25))) + 
  xlab("Bureaucratic quality") + ylab("Effect of Treaty institutionalization")

ggsave(plot = ols.bq.abgmm.mep, filename = "Output/megg_ols_bq_abgmm.pdf", width = 4, height = 4)

instcoop <- seq(0, 4, by = 0.2)
b1.bq <- ols.abgmm$coefficients[2]
b3.interact <- ols.abgmm$coefficients[9]
var1.bq <- plm::vcovHC(ols.abgmm)[2, 2]
var3.interact <- plm::vcovHC(ols.abgmm)[9, 9]
cov13.bq_interact <- plm::vcovHC(ols.abgmm)[2, 9]
mf.bq <- b1.bq + b3.interact * instcoop
var.bq <- var1.bq + instcoop^2 * var3.interact + 2 * instcoop * cov13.bq_interact
se.bq <- sqrt(var.bq)
low1.bq <- mf.bq - 1.96 * se.bq
up1.bq <- mf.bq + 1.96 * se.bq

ols.abgmm.plot.dat <- data.frame(s2 = instcoop, eff1 = mf.bq, se.eff1 = se.bq, low1 = low1.bq, up1 = up1.bq)

ols.instcoop.abgmm.mep <- ggplot(data = ols.abgmm.plot.dat, aes(x = s2, y = eff1)) + 
  geom_hline(yintercept = 0, linetype = "dashed") + theme_jk() +
  geom_ribbon(aes(x = s2, ymin = low1, ymax = up1), alpha = 0.25, color = NA) +
  geom_line() + 
  scale_x_continuous(breaks = c(seq(0, 4, by = 1))) + 
  xlab("Treaty institutionalization") + ylab("Effect of Bureaucratic Quality")

ggsave(plot = ols.instcoop.abgmm.mep, filename = "Output/megg_ols_inst_abgmm.pdf", width = 4, height = 4)

###  7a. Year FE
# ols.bq.yfe

ols.bq.yfe.mep <- ggintfun(ols.bq.yfe, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE) + scale_y_continuous(limits = c(-0.05, 0.25))
ggsave(plot = ols.bq.yfe.mep, filename = "Output/megg_ols_bq_yfe.pdf", width = 4, height = 4)

ols.inst.yfe.mep <- ggintfun(ols.bq.yfe, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE)
ggsave(plot = ols.inst.yfe.mep, filename = "Output/megg_ols_inst_yfe.pdf", width = 4, height = 4)


###  8. Dyad & year RE
# ols.bq.mlm

ols.bq.mlm.mep <- ggintfun(ols.bq.mlm, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE) + scale_y_continuous(limits = c(-0.05, 0.25))
ggsave(plot = ols.bq.mlm.mep, filename = "Output/megg_ols_bq_mlm.pdf", width = 4, height = 4)

ols.inst.mlm.mep <- ggintfun(ols.bq.mlm, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE)
ggsave(plot = ols.inst.mlm.mep, filename = "Output/megg_ols_inst_mlm.pdf", width = 4, height = 4)


###  9. Rivalry
# ols.bq.riv

ols.bq.riv.mep <- ggintfun(ols.bq.riv, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE) + scale_y_continuous(limits = c(-0.05, 0.25))
ggsave(plot = ols.bq.riv.mep, filename = "Output/megg_ols_bq_riv.pdf", width = 4, height = 4)

ols.inst.riv.mep <- ggintfun(ols.bq.riv, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE)
ggsave(plot = ols.inst.riv.mep, filename = "Output/megg_ols_inst_riv.pdf", width = 4, height = 4)


###  10. IV
# obtain coefficient estimates from Stata output

# Instrumental variables (2SLS) regression               Number of obs =   10981
# Wald chi2(9)  =  215.27
# Prob > chi2   =  0.0000
# R-squared     =  0.0202
# Root MSE      =  .84477
# 
# --------------------------------------------------------------------------------------
#             bar_avg2 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
# ---------------------+----------------------------------------------------------------
#        instcoop_max2 |   .0624967   .0134649     4.64   0.000     .0361059    .0888874
#     instcoopXburqual |   .1137409   .0234725     4.85   0.000     .0677357    .1597461
# icrg_burqual_min_std |   .0504825   .0411605     1.23   0.220    -.0301907    .1311557
#         wres_min_log |  -.0049844   .0054853    -0.91   0.364    -.0157354    .0057665
#           ntreaties2 |  -.0150114   .0049997    -3.00   0.003    -.0248107   -.0052122
#             dyaddem7 |  -.0888734   .0252874    -3.51   0.000    -.1384358   -.0393109
#        gdpmaxPWT_log |  -.0051645    .008524    -0.61   0.545    -.0218712    .0115422
#         lnpowerratio |   .0126569   .0056653     2.23   0.025      .001553    .0237607
#             alliance |  -.0016987   .0195147    -0.09   0.931    -.0399469    .0365494
#                _cons |   .1890805   .0793796     2.38   0.017     .0334994    .3446616
# --------------------------------------------------------------------------------------
# Instrumented:  instcoop_max2 instcoopXburqual
# Instruments:   icrg_burqual_min_std wres_min_log ntreaties2 dyaddem7 gdpmaxPWT_log
#                lnpowerratio alliance instcoop_max2_IV instcoopXburqual_IV

# . matrix list e(V)
# 
# symmetric e(V)[10,10]
#               instcoop_m~2  instcoopXb~l  icrg_burqu~d  wres_min_log    ntreaties2      dyaddem7  gdpmaxPWT_~g  lnpowerratio      alliance         _cons
# instcoop_m~2      .0001813
# instcoopXb~l     -.0002064     .00055096
# icrg_burqu~d     .00016032    -.00043829     .00169419
# wres_min_log     9.998e-06    -3.151e-06     .00003211     .00003009
#   ntreaties2    -.00002837    -.00001894     8.925e-07    -3.076e-06       .000025
#     dyaddem7    -4.715e-06    -.00002155    -.00024025    -.00002402    -3.259e-06     .00063945
# gdpmaxPWT_~g     3.386e-06    -.00001064    -.00011032    -5.376e-07     2.450e-06    -.00008722     .00007266
# lnpowerratio     4.442e-06     .00001368    -.00004317    -1.813e-06    -3.008e-06    -2.406e-06    -2.980e-06      .0000321
#     alliance    -9.024e-06     .00001786     .00006026    -.00001334    -.00001333    -.00004412     8.344e-06     8.360e-06     .00038083
#        _cons    -.00020488     .00025833     .00015497    -.00022912     9.411e-06     .00090446    -.00054888    -4.797e-06    -.00007665     .00630112

icrg_burqual_min_std <- seq(0, 1, by = 0.05)
b1.instcoop <- .0624967
b3.interact <- .1137409
var1.instcoop <- .0001813
var3.interact <- .00055096
cov13.instcoop_interact <- -.0002064
mf.instcoop <- b1.instcoop + b3.interact * icrg_burqual_min_std
var.instcoop <- var1.instcoop + icrg_burqual_min_std^2 * var3.interact + 2 * icrg_burqual_min_std * cov13.instcoop_interact
se.instcoop <- sqrt(var.instcoop)
low1.instcoop <- mf.instcoop - 1.96 * se.instcoop
up1.instcoop <- mf.instcoop + 1.96 * se.instcoop

iv.plot.dat <- data.frame(s2 = icrg_burqual_min_std, eff1 = mf.instcoop, se.eff1 = se.instcoop, low1 = low1.instcoop, up1 = up1.instcoop)

ols.bq.iv.mep <- ggplot(data = iv.plot.dat, aes(x = s2, y = eff1)) + 
          geom_hline(yintercept = 0, linetype = "dashed") + theme_jk() +
          geom_ribbon(aes(x = s2, ymin = low1, ymax = up1), alpha = 0.25, color = NA) +
          geom_line() + 
          scale_x_continuous(breaks = c(seq(0, 1, by = 0.25))) + 
          scale_y_continuous(limits = c(-0.05, 0.25)) +
          xlab("Bureaucratic quality") + ylab("Effect of Treaty institutionalization")

ggsave(plot = ols.bq.iv.mep, filename = "Output/megg_ols_bq_iv.pdf", width = 4, height = 4)

###  11. Breaking up institutions

# ols.bq.institution
# ols.bq.monitoring
# ols.bq.enforcement
# ols.bq.conflictres


b1.inst.inst <- coef(ols.bq.institution)[2]
b1.moni.moni <- coef(ols.bq.monitoring)[2]
b1.enfo.enfo <- coef(ols.bq.enforcement)[2]
b1.conf.conf <- coef(ols.bq.conflictres)[2]

b3.inst.inst_bq <- coef(ols.bq.institution)[10]
b3.moni.moni_bq <- coef(ols.bq.monitoring)[10]
b3.enfo.enfo_bq <- coef(ols.bq.enforcement)[10]
b3.conf.conf_bq <- coef(ols.bq.conflictres)[10]

var.b1.inst.inst <- vcov(ols.bq.institution)[2, 2]
var.b1.moni.moni <- vcov(ols.bq.monitoring)[2, 2]
var.b1.enfo.enfo <- vcov(ols.bq.enforcement)[2, 2]
var.b1.conf.conf <- vcov(ols.bq.conflictres)[2, 2]

var.b3.inst.inst_bq <- vcov(ols.bq.institution)[10, 10]
var.b3.moni.moni_bq <- vcov(ols.bq.monitoring)[10, 10]
var.b3.enfo.enfo_bq <- vcov(ols.bq.enforcement)[10, 10]
var.b3.conf.conf_bq <- vcov(ols.bq.conflictres)[10, 10]

cov.b1b3.inst <- vcov(ols.bq.institution)[2, 10]
cov.b1b3.moni <- vcov(ols.bq.monitoring)[2, 10]
cov.b1b3.enfo <- vcov(ols.bq.enforcement)[2, 10]
cov.b1b3.conf <- vcov(ols.bq.conflictres)[2, 10]

# Simulate data for different values of moderator

sim_icrg_burqual_min_std <- seq(0, 1, by = 0.1)

# Calculate conditional effects and their SE

eff.inst <- b1.inst.inst + b3.inst.inst_bq * sim_icrg_burqual_min_std
eff.moni <- b1.moni.moni + b3.moni.moni_bq * sim_icrg_burqual_min_std
eff.enfo <- b1.enfo.enfo + b3.enfo.enfo_bq * sim_icrg_burqual_min_std
eff.conf <- b1.conf.conf + b3.conf.conf_bq * sim_icrg_burqual_min_std

eff.inst.se <- sqrt(var.b1.inst.inst + sim_icrg_burqual_min_std^2 * var.b3.inst.inst_bq + 2 * sim_icrg_burqual_min_std * cov.b1b3.inst)
eff.moni.se <- sqrt(var.b1.moni.moni + sim_icrg_burqual_min_std^2 * var.b3.moni.moni_bq + 2 * sim_icrg_burqual_min_std * cov.b1b3.moni)
eff.enfo.se <- sqrt(var.b1.enfo.enfo + sim_icrg_burqual_min_std^2 * var.b3.enfo.enfo_bq + 2 * sim_icrg_burqual_min_std * cov.b1b3.enfo)
eff.conf.se <- sqrt(var.b1.conf.conf + sim_icrg_burqual_min_std^2 * var.b3.conf.conf_bq + 2 * sim_icrg_burqual_min_std * cov.b1b3.conf)

eff.inst.lo <- eff.inst - 1.96 * eff.inst.se
eff.moni.lo <- eff.moni - 1.96 * eff.moni.se
eff.enfo.lo <- eff.enfo - 1.96 * eff.enfo.se
eff.conf.lo <- eff.conf - 1.96 * eff.conf.se

eff.inst.hi <- eff.inst + 1.96 * eff.inst.se
eff.moni.hi <- eff.moni + 1.96 * eff.moni.se
eff.enfo.hi <- eff.enfo + 1.96 * eff.enfo.se
eff.conf.hi <- eff.conf + 1.96 * eff.conf.se

# Plots

inst_sep_d <- data.frame(sim_icrg_burqual_min_std = rep(sim_icrg_burqual_min_std, 4),
                         effect = c(eff.inst, eff.moni, eff.enfo, eff.conf),
                         lo = c(eff.inst.lo, eff.moni.lo, eff.enfo.lo, eff.conf.lo),
                         hi = c(eff.inst.hi, eff.moni.hi, eff.enfo.hi, eff.conf.hi),
                         inst = c(rep("Treaty delegates to IGO", length(sim_icrg_burqual_min_std)), rep("Treaty provides Monitoring", length(sim_icrg_burqual_min_std)), rep("Treaty provides Enforcement", length(sim_icrg_burqual_min_std)), rep("Treaty provides Conflict Resolution", length(sim_icrg_burqual_min_std))))

me_inst_sep_all <- ggplot(data = inst_sep_d, aes(x = sim_icrg_burqual_min_std, y = effect)) + 
  geom_hline(yintercept = 0, linetype = "dashed") + theme_jk() +
  geom_ribbon(aes(x = sim_icrg_burqual_min_std, ymin = lo, ymax = hi), alpha = 0.25, color = NA) +
  geom_line() + 
  scale_x_continuous(breaks = c(seq(0, 1, by = 0.25))) + 
  scale_y_continuous(limits = c(-0.4, 0.55)) +
  xlab("Bureaucratic quality") + ylab("Effect of Treaty institutionalization") +
  facet_wrap( ~ inst, scales = "free_y")

ggsave(plot = me_inst_sep_all, filename = "Output/megg_inst_sep_all.pdf", width = 8, height = 8)

mod <- c("Treaty delegates to IGO", "Treaty provides Monitoring", "Treaty provides Enforcement", "Treaty provides Conflict Resolution")
modshort <- c("IGO", "Monitoring", "Enforcement", "ConflictResolution")
for(i in 1:length(mod)){
  p <- ggplot(data = inst_sep_d[inst_sep_d$inst == paste(mod[i]), ], aes(x = sim_icrg_burqual_min_std, y = effect)) + 
    geom_hline(yintercept = 0, linetype = "dashed") + theme_jk() +
    geom_ribbon(aes(x = sim_icrg_burqual_min_std, ymin = lo, ymax = hi), alpha = 0.25, color = NA) +
    geom_line() + 
    scale_x_continuous(breaks = c(seq(0, 1, by = 0.25))) + 
    scale_y_continuous(limits = c(-0.4, 0.55)) +
    xlab("Bureaucratic quality") + ylab(paste(mod[i],":\nConditional effect", sep = ""))
  
  ggsave(plot = p, filename = paste("Output/megg_inst_sep_",modshort[i],".pdf", sep = ""), width = 4, height = 4)
}


###  13. Control for shared IGO memberships
# ols.bq.igos

ols.bq.igos.mep <- ggintfun(ols.bq.igos, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE) + scale_y_continuous(limits = c(-0.05, 0.25))
ggsave(plot = ols.bq.igos.mep, filename = "Output/megg_ols_bq_igos.pdf", width = 4, height = 4)

ols.inst.igos.mep <- ggintfun(ols.bq.igos, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE)
ggsave(plot = ols.inst.igos.mep, filename = "Output/megg_ols_inst_igos.pdf", width = 4, height = 4)


###  15. States in treaty
# ols.bq.tm

ols.bq.tm.mep <- ggintfun(ols.bq.tm, varnames = c("instcoop_max2", "icrg_burqual_min_std"), varlabs = c("Treaty institutionalization", "Bureaucratic quality"), rug = FALSE) + scale_y_continuous(limits = c(-0.05, 0.25))
ggsave(plot = ols.bq.tm.mep, filename = "Output/megg_ols_bq_tm.pdf", width = 4, height = 4)

ols.inst.tm.mep <- ggintfun(ols.bq.tm, varnames = c("icrg_burqual_min_std", "instcoop_max2"), varlabs = c("Bureaucratic quality", "Treaty institutionalization"), rug = FALSE)
ggsave(plot = ols.inst.tm.mep, filename = "Output/megg_ols_inst_tm.pdf", width = 4, height = 4)



######################################
### Summary statistics
######################################

# Table

d_sum <- dat[ , c("bar_avg2", "bar_avg", "bar_pos2only", "instcoop_max2", "instcoop_max_05", "institution_max2", "monitoring_max2", "enforcement_max2", "conflictres_max2", "icrg_burqual_min_std", "wres_min_log", "ntreaties2", "dyaddem7", "gdpmaxPWT_log", "lnpowerratio", "alliance", "rivalry_th", "sharedIGOcount", "treatymembers.max2", "year")]
d_sum$sample <- ifelse(is.na(d_sum$bar_avg2) == FALSE & is.na(d_sum$instcoop_max2) == FALSE & is.na(d_sum$icrg_burqual_min_std) == FALSE & is.na(d_sum$wres_min_log) == FALSE & is.na(d_sum$ntreaties2) == FALSE & is.na(d_sum$dyaddem7) == FALSE & is.na(d_sum$gdpmaxPWT_log) == FALSE & is.na(d_sum$lnpowerratio) == FALSE & is.na(d_sum$alliance) == FALSE, 1, 0)
d_sum <- d_sum[d_sum$sample == 1, ]
d_sum$sample <- NULL
d_sum$year <- NULL

Mean <- apply(X = d_sum, MARGIN = 2, FUN = function(x) mean(x, na.rm = TRUE))
SD <- apply(X = d_sum, MARGIN = 2, FUN = function(x) sd(x, na.rm = TRUE))
Min <- apply(X = d_sum, MARGIN = 2, FUN = function(x) min(x, na.rm = TRUE))
Max <- apply(X = d_sum, MARGIN = 2, FUN = function(x) max(x, na.rm = TRUE))
N <- apply(X = d_sum, MARGIN = 2, FUN = function(x) length(x[is.na(x) == FALSE]))

Variable <- c("Water-related cooperation (all years)", "Water-related cooperation (years with events only)", "At least one cooperative event", "Treaty institutionalization", "Treaty institutionalization (no treaties separate)", "Treaty created IGO", "Treaty provides Monitoring", "Treaty provides Enforcement", "Treaty provides conflict resolution", "Bureaucratic quality (lower)", "Water availability (lower)", "Treaty count", "Democratic dyad", "GDP p.c. (higher, logged)", "Power ratio", "Alliance", "Rivalry", "Shared IGO memberships", "Number of treaty members")

sumtab <- data.frame(Variable, Mean, SD, Min, Max, N)
sumtab_digits <- matrix(c(rep(2, nrow(sumtab)), rep(2, nrow(sumtab)), rep(2, nrow(sumtab)), rep(2, nrow(sumtab)), ifelse(sumtab$Min == 0, 0, 2), ifelse(sumtab$Min == 0, 0, 2), rep(2, nrow(sumtab))), nrow = nrow(sumtab), ncol = ncol(sumtab) + 1, byrow = FALSE)
sumtab_x <- xtable(sumtab, caption = c("Summary statistics"), label = "tab:sumstats", align = "llccccc", digits = sumtab_digits)

print(sumtab_x, type = "latex", file = "Output/tab_sumstats.tex", include.rownames = FALSE, include.colnames = TRUE, caption.placement = "top", size = "scriptsize", hline.after = c(-1,0,nrow(sumtab_x)), booktabs = TRUE)


######################################
### Maps of River treaty counts (Figure A2)
######################################

# Run this first before the rest

library("ggplot2")
library("maptools")
reg <- map_data("world")
P1 <- ggplot(reg, aes(long, lat)) + 
  geom_polygon(aes(group = group), show.legend = FALSE)
worldmap <- map("world", fill = TRUE, col = "transparent", plot = FALSE)
IDs <- sapply(strsplit(worldmap$names, ":"), function(x) x[1])
world_sp <- map2SpatialPolygons(worldmap, IDs=IDs, proj4string=CRS("+proj=longlat +datum=wgs84"))

library("rgdal")
library("countrycode")
library("foreign")

gpclibPermit()

# Load shapefile

world.map <- readOGR(dsn = "TM_WORLD_BORDERS_SIMPL-0.3", layer = "TM_WORLD_BORDERS_SIMPL-0.3")

# Remove Antarctica
world.map <- world.map[world.map$NAME != "Antarctica", ]

# Prepare for ggplot2
world.ggmap <- fortify(world.map, region = "NAME")

# Just to be sure
world.ggmap <- world.ggmap[world.ggmap$id != "Antarctica", ]

# Supply COW codes
world.ggmap$cowcode <- countrycode(world.ggmap$id, origin = "country.name", destination = "cown")

# See if any COW codes are missing
sort(unique(world.ggmap[is.na(world.ggmap$cowcode)==TRUE, ]$id) )
world.ggmap[world.ggmap$id == "Korea, Democratic People's Republic of", ]$cowcode <- 731
world.ggmap[world.ggmap$id == "Serbia", ]$cowcode <- 345

# Read in RT data by country
rt.dat <- read.csv("rtd_treaty-country-year.csv")
rt.dat <- rt.dat[, c("cowcode", "year", "treaty_metadata_docid", "treatyinst.ind", "treatyinst.ind2")]
rt.dat$treaty <- c(rep(1, length(rt.dat$cowcode)))
library(dplyr)
rt.dat.country <- summarise(group_by(rt.dat, cowcode, year), treatycount = sum(treaty), treatyinst.ind.max = max(treatyinst.ind), treatyinst.ind2.max = max(treatyinst.ind2), treatyinst.ind.min = min(treatyinst.ind), treatyinst.ind2.min = min(treatyinst.ind2))
rt.dat.country.2000 <- rt.dat.country[rt.dat.country$year == 2000, ]
rt.dat.country.2000 <- as.data.frame(rt.dat.country.2000)
rt.dat.country.2000$cowcode <- as.numeric(rt.dat.country.2000$cowcode)

# Use plyr to merge
map.data <- left_join(world.ggmap, rt.dat.country.2000, by = "cowcode")

# Center points of countries for later
centroids <- summarise(group_by(map.data, id), long.center = mean(long), lat.center = mean(lat))
# Note: a better solution would be capitals
# also try the coordinates() command, via ugly workaround

# Data frame for plotting
df <- data.frame(id = unique(map.data$id), cowcode = tapply(map.data$cowcode, map.data$id, mean), treatycount = tapply(map.data$treatycount, map.data$id, mean), treatyinst.ind2.max = tapply(map.data$treatyinst.ind2.max, map.data$id, mean), treatyinst.ind.max = tapply(map.data$treatyinst.ind.max, map.data$id, mean))
# df <- data.frame(df, centroids.df)
# df <- inner_join(df[is.na(df$cowcode) == FALSE, ], centroids.df, by = "cowcode")

# Delete countries w/out COW codes
df <- df[is.na(df$cowcode) == FALSE, ]

# Need to think better about color choices. RColorBrewer
library("RColorBrewer")

p2000 <- ggplot() + geom_map(aes(fill = treatycount, map_id = id), data = df, map = world.ggmap, colour = "black", size = 0.125) + expand_limits(x = world.ggmap$long, y = world.ggmap$lat) + scale_fill_gradient(low = "gray90", high = "gray10", na.value = "lightyellow", name = "Count of \nriver treaties") + coord_equal() + theme_bw() + xlab("") + ylab("") +  theme(line = element_blank(), axis.text = element_blank(), legend.position = c(0.1, 0.35))

pdf("Output/rt2000_map.pdf", width = 11, height = 11)
p2000
dev.off()

######################################
### End of file
######################################

sessionInfo()

# R version 3.4.2 (2017-09-28)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS Sierra 10.12.6
# 
# Matrix products: default
# BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] RColorBrewer_1.1-2      rgdal_1.2-15            maps_3.2.0              maptools_0.9-2         
# [5] sp_1.2-5                systemfit_1.1-20        lmtest_0.9-35           zoo_1.8-0              
# [9] car_2.1-5               RStata_1.1.1            Formula_1.2-2           directlabels_2017.03.31
# [13] bindrcpp_0.2            countrycode_0.19        rio_0.5.5               ggplot2_2.2.1          
# [17] effects_4.0-0           carData_3.0-0           gridExtra_2.3           dplyr_0.7.4            
# [21] xtable_1.8-2            lme4_1.1-14             Matrix_1.2-11           texreg_1.36.23         
# [25] foreign_0.8-69         
# 
# loaded via a namespace (and not attached):
#   [1] Rcpp_0.12.13        bdsmatrix_1.3-2     lattice_0.20-35     assertthat_0.2.0   
# [5] digest_0.6.12       R6_2.2.2            cellranger_1.1.0    plyr_1.8.4         
# [9] MatrixModels_0.4-1  survey_3.32-1       rlang_0.1.2         lazyeval_0.2.1     
# [13] curl_3.0            readxl_1.0.0        minqa_1.2.4         data.table_1.10.4-3
# [17] SparseM_1.77        nloptr_1.0.4        labeling_0.3        splines_3.4.2      
# [21] readr_1.1.1         bit_1.1-12          munsell_0.4.3       compiler_3.4.2     
# [25] pkgconfig_2.0.1     mgcv_1.8-22         nnet_7.3-12         tibble_1.3.4       
# [29] quadprog_1.5-5      MASS_7.3-47         grid_3.4.2          nlme_3.1-131       
# [33] gtable_0.2.0        magrittr_1.5        scales_0.5.0        plm_1.6-5          
# [37] sandwich_2.4-0      openxlsx_4.0.17     tools_3.4.2         forcats_0.2.0      
# [41] bit64_0.9-7         glue_1.2.0          hms_0.3             pbkrtest_0.4-7     
# [45] parallel_3.4.2      survival_2.41-3     yaml_2.1.14         colorspace_1.3-2   
# [49] gpclib_1.5-5        bindr_0.1           haven_1.1.0         quantreg_5.34